{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GO enrichment for retained paralogs example\n",
    "\n",
    "Here I describe how to get to a GO enrichment analysis for retained paralogs from a WGD event. To this end I use the `wgd` library to perform mixture modeling and select paralogs from a particular mixture component. The GO enrichment is then performed using the `goenrich` library (https://github.com/jdrudolph/goenrich). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from wgd.viz import plot_selection\n",
    "from wgd.modeling import fit_bgmm, get_component_probabilities, filter_group_data, plot_all_models_bgmm, get_array_for_mixture\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import goenrich"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Get the distributions\n",
    "\n",
    "The distributions were computed using the '`wgd ks`' and '`wgd syn`' commands from the '`wgd`' CLI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDAAAAKUCAYAAADyycjUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xu4ZHV95/v3h3u8ILeIDN0RDB0TFEW6B83RyaCoEaJCTlAxibZI0jrRJMYkiuE5o4m56CQT4+3gECGCURFRh47BKEGNxxxRuwUaBJW21dAt0qKAooIi3/mjfo1Fu/fufaldtarq/Xqe9dRav3XZ31W1d+1ffet3SVUhSZIkSZLUZbuNOgBJkiRJkqRdMYEhSZIkSZI6zwSGJEmSJEnqPBMYkiRJkiSp80xgSJIkSZKkzjOBIUmSJEmSOs8EhiRJkiRJ6jwTGJIkSZIkqfNMYEiSJEmSpM7bY9QBLLMadQCSJE2gjDqAEbJuIUnS4M2rbjH0FhhJHprkyr7l20lekuSAJJcmub497t+OT5I3JNmcZFOSY4YdsyRJ6lm5ciWL/V+c5IYkj2rrSfLVJEcMNkJJkjSphp7AqKovVNXRVXU0sBr4HvB+4AzgsqpaBVzWtgFOAFa1ZR1w1rBjliRp3N111138wi/8AqtWrVrQeStWrOCKK64A4JZbbuFrX/sawHUL/fnti4kHAdcCVM+Dq2rzQq8lSZKm06jHwDge+FJVfRU4CTivlZ8HnNzWTwLObxWdy4H9khwy/FAlSRpfb3nLW9i+fTtbtmzhu9/97rzOufnmm7nppps48sgjAbj66qs57LDDqKrvLyKEo4AvVtWdizhXkiRp5AmMU4F3tfWDq+rGtv514OC2fihwQ985W1uZJEmah9tuu40//dM/5U1vehO7774711xzzS7P2bx5MytXruTuu+/mwAMP5MADD2TTpk0cccQRtK6d30jytSRP2nFOkgcmWZ/kpiTfSfJPSfZtux8BXN137OlJ/veuzkvygiSXJHlzkpsX+DMlSdIEGVkCI8lewNOB9+y8r6qKBQ6SlWRdkg1t2XXNTJKkKfHqV7+ahzzkITz72c/m537u59i0adMuzzniiCP4m7/5G0455RRuv/12vvnNb3L11VezYcMGgH+m90XD/wJe3nfavsAbgZ8BHgwcBLyg7TsK6P/BjwSumsd5jwQeA6wHHrjAnylJkibIKFtgnAB8tqpuats37ega0h63t/JtwMq+81a0snupqrOrak1VrQHuWL6wJUkaH1u2bOGNb3wjr33tawF42MMeNq8EBsBVV13F0Ucffc/2pk2bOPPMM6mqD1XV3bTxLHaoqs1VdWlV3VlV3wIuBfZvu+/VAoO+BMY8znvNIn+mJEmaIKNMYDybH3cfgd43K2vb+lrg4r7y57bRyh8D3NbX1USSJM3hZS97GU984hM57rjjgF4C46qrrpr7pObKK6/kkY98JABVxTXXXMPTnva0/kMeTl9CIckzkvx7ku1JbqU3IPcXk6Qd2585eQQtgbGL844C/mmhP3NeNyhJksbKHqP4oUnuCzyJezfxfA1wYZLTga8Cz2zllwAnApvpzVhy2hBDlSRpbH3iE5/gve99L/vuuy8PetCDAPj+97/Pbrvt+vuLu+++m2uuueaeFhhf/vKXgV7Xkj6PAnaMY/EE4LXAs4Ar2v6vAFcChwM/aoN2k+TB9OogW3Zx3mHtuC8s4mdKkqQJM5IWGFX13ao6sKpu6yv7ZlUdX1WrquqJrRnojmnWXlRVP1tVR1XVhlHELEnSOKkqXvrSl/LCF76QL3zhC1x55ZVceeWVfPjDH+bWW2/lP/7jP+Y8//vf/z7f//73ufvuu4Fe95GjjjqKXqOIezyKH49j8Uh6g25fRa8Lx7n0xqy4ll5ri/7xqR4JbGpjXu3qvKtb15GF/kxJkjRhRj0LiSRJWgbveMc7uOmmm/jrv/5rHvSgB92zHHvssdz//ve/ZxyME044gb/8y7/8ifPve9/78sIXvpAjjzySFStWcPXVV9/TnQQgyUHAg/hxYuIdwJ7At4APANcD11bVD5h5AM9N8zjvEfS1pljgz5QkSRMmvS8/JkuSDW0wz8m7OUmSRi+7PmRiWbeQJGnw5lW3sAWGJEmSJEnqPBMYkiRp7CXZL8lFST6f5Lokv5jkgCSXJrm+Pe7fjk2SNyTZnGRTkmNGHb8kSdq1kcxCom7YuHHjrPtWr149xEgkSVqy1wP/UlWnJNkLuA/wJ8BlVfWaJGfQm2L15cAJwKq2PBo4qz1KkjTWJv0zngmMCTTXLy1Mxi+uJEk7JHkA8EvA8wDaIJ4/SHIScFw77DzgY/QSGCcB57dZUC5vrTcOqaobhxy6JElaALuQSJKkcXc48A3gH5JckeStSe4LHNyXlPg6cHBbP5Te9Ks7bG1lkiSpw0xgSJKkcbcHcAxwVlU9Cvguve4i92itLRY0g0iSdUk2tOWaXZ8hSZKWkwkMSZI07rYCW6vqU237InoJjZuSHALQHre3/duAlX3nr2hl91JVZ1fVmjY1+x3LFbwkSZofExiSJGmsVdXXgRuSPLQVHQ9cC6wH1raytcDFbX098Nw2G8ljgNsc/0KSpO5zEE9JkjQJfhd4R5uBZAtwGr0vai5McjrwVeCZ7dhLgBOBzcD32rGSJKnjTGBIkqSxV1VXAmtm2HX8DMcW8KJlD0qSJA2UXUgkSZIkSVLnmcCQJEmSJEmdZwJDkiRJkiR1ngkMSZIkSZLUeSYwJEmSJElS5zkLiSRJkiRJU2Ljxo1z7l+9evWQIlk4W2BIkiRJkqTOM4EhSZIkSZI6zy4kY2au5j5dbuojSZIkSdJS2AJDkiRJkiR1ngkMSZIkSZLUeSYwJEmSJElS55nAkCRJkiRJnWcCQ5IkSZIkdZ4JDEmSJEmS1HkmMCRJkiRJUueZwJAkSZIkSZ1nAkOSJEmSJHXeSBIYSfZLclGSzye5LskvJjkgyaVJrm+P+7djk+QNSTYn2ZTkmFHELEmSJEmSRmdULTBeD/xLVf088EjgOuAM4LKqWgVc1rYBTgBWtWUdcNbww5UkSZIkSaM09ARGkgcAvwScA1BVP6iqW4GTgPPaYecBJ7f1k4Dzq+dyYL8khww5bEmSJEmSNEKjaIFxOPAN4B+SXJHkrUnuCxxcVTe2Y74OHNzWDwVu6Dt/ayuTJEmSJElTYhQJjD2AY4CzqupRwHf5cXcRAKqqgFrIRZOsS7IhyQbgoEEFK0mSJEmSRm8UCYytwNaq+lTbvoheQuOmHV1D2uP2tn8bsLLv/BWt7F6q6uyqWlNVa4Cblyt4SZIkSZI0fENPYFTV14Ebkjy0FR0PXAusB9a2srXAxW19PfDcNhvJY4Db+rqaSJIkSZKkKbDHiH7u7wLvSLIXsAU4jV4y5cIkpwNfBZ7Zjr0EOBHYDHyvHash2bhx46z7Vq9ePcRIJEmSJEnTbCQJjKq6Elgzw67jZzi2gBcte1CSJEmSJKmzRtUCQ5IkSZI0B1tDS/dmAkOSJEmSpA6bK5kF05PQGsUsJJIkSZIkSQtiAkOSJEmSJHWeCQxJkiRJktR5joEhSZIkSdKIOL7F/NkCQ5IkSZIkdZ4JDEmSNPaSfCXJ1UmuTLKhlR2Q5NIk17fH/Vt5krwhyeYkm5IcM9roJUnSfJjAkCRJk+LxVXV0Va1p22cAl1XVKuCytg1wArCqLeuAs4YeqSRJWjATGJIkaVKdBJzX1s8DTu4rP796Lgf2S3LIKAKUJEnzZwJDkiRNggI+nGRjknWt7OCqurGtfx04uK0fCtzQd+7WViZJkjrMWUgkSdIkeFxVbUvyQODSJJ/v31lVlaQWcsGWCNmRDDloQHFKkqRFsgWGJEkae1W1rT1uB94PHAvctKNrSHvc3g7fBqzsO31FK9v5mmdX1Zo2psbNyxi+JEmaBxMYkiRprCW5b5L771gHngxcA6wH1rbD1gIXt/X1wHPbbCSPAW7r62oiSZI6yi4kkiRp3B0MvD8J9Oo276yqf0nyGeDCJKcDXwWe2Y6/BDgR2Ax8Dzht+CFLkqSFMoEhSZLGWlVtAR45Q/k3geNnKC/gRUMITZIkDZBdSCRJkiRJUueZwJAkSZIkSZ1nAkOSJEmSJHWeY2BIkiRJ0hBt3Lhxzv2rV68eUiTSeLEFhiRJkiRJ6jwTGJIkSZIkqfNMYEiSJEmSpM5zDAxJkiRJmmCOuaFJYQsMSZIkSZLUeSYwJEmSJElS59mFRJIkSZLGlN1DNE1MYGjJfNOUJEmSJC03u5BIkiRJkqTOM4EhSZIkSZI6byRdSJJ8BfgO8CPgrqpak+QA4N3AYcBXgGdW1S1JArweOBH4HvC8qvrsKOKWJEmSJGm+5upub1f7hRtlC4zHV9XRVbWmbZ8BXFZVq4DL2jbACcCqtqwDzhp6pJIkSZIkaaS61IXkJOC8tn4ecHJf+fnVczmwX5JDRhGgJEmSJEkajVElMAr4cJKNSda1soOr6sa2/nXg4LZ+KHBD37lbW5kkSZIkSZoSo5pG9XFVtS3JA4FLk3y+f2dVVZJayAVbImRHMuSgAcUpSZIkSZI6YCQtMKpqW3vcDrwfOBa4aUfXkPa4vR2+DVjZd/qKVrbzNc+uqjVtTI2blzF8SZIkSZI0ZENvgZHkvsBuVfWdtv5k4M+A9cBa4DXt8eJ2ynrgxUkuAB4N3NbX1WTkHFVWkiRJkqTlN4ouJAcD7+/NjsoewDur6l+SfAa4MMnpwFeBZ7bjL6E3hepmetOonjb8kCVJkiRJ0igNPYFRVVuAR85Q/k3g+BnKC3jREEKTJEmSJEkd1aVpVCVJkiRJkmZkAkOSJEmSJHWeCQxJkiRJktR5oxjEU5IkSZKksTXXbJTgjJTLxRYYkiRJkiSp80xgSJIkSZKkzjOBIUmSJEmSOs8xMDpkrn5U9qGSJEmSJE0zW2BIkiRJkqTOswWGJEmSJEm6ly72EDCBIUmSJEkD0sUPfdKksAuJJEmSJEnqPBMYkiRpIiTZPckVST7Qtg9P8qkkm5O8O8lerXzvtr257T9slHFLkqT5sQuJJEmaFL8PXAfs27ZfC7yuqi5I8hbgdOCs9nhLVR2R5NR23LNGEbCk8TFX1xCwe4g0DLbAkCRJYy/JCuBXgLe27QBPAC5qh5wHnNzWT2rbtP3Ht+MlSVKHmcCQJEmT4O+AlwF3t+0DgVur6q62vRU4tK0fCtwA0Pbf1o6XJEkdZgJDkiSNtSRPBbZX1dztuxd+3XVJNiTZABw0yGtLkqSFcwwMSZI07h4LPD3JicA+9MbAeD2wX5I9WiuLFcC2dvw2YCWwNckewAOAb+580ao6GzgboCUxJEnSCJnAGAIH/JEkaflU1SuAVwAkOQ74o6r6jSTvAU4BLgDWAhe3U9a37U+2/R+pqhp23JIkaWHsQiJJkibVy4GXJtlMb4yLc1r5OcCBrfylwBkjik+SJC2ALTAkSdLEqKqPAR9r61uAY2c45g7gGUMNTJIkLZktMCRJkiRJUueZwJAkSZIkSZ1nAkOSJEmSJHWeCQxJkiRJktR5JjAkSZIkSVLnmcCQJEmSJEmdt+QERpKfTbJ3Wz8uye8l2W/poUmSpGljvUKSJM1mEC0w3gv8KMkRwNnASuCdA7iuJEmaPtYrJEnSjAaRwLi7qu4CfhV4Y1X9MXDIAK4rSZKmj/UKSZI0o0EkMH6Y5NnAWuADrWzPXZ2UZPckVyT5QNs+PMmnkmxO8u4ke7Xyvdv25rb/sAHELEmSumlR9QpJkjT5BpHAOA34ReAvqurLSQ4H3j6P834fuK5v+7XA66rqCOAW4PRWfjpwSyt/XTtOkiRNpsXWKyRJ0oRbcgKjqq4FXg58tm1/uarmTDIkWQH8CvDWth3gCcBF7ZDzgJPb+kltm7b/+Ha8JEmaMIupV0iSpOmwx1IvkORpwN8AewGHJzka+LOqevocp/0d8DLg/m37QODW1ucVYCtwaFs/FLgBoKruSnJbO/7mpcYuSZK6ZZH1CknSAGzcuHHWfatXrx5iJKPl89Bdg+hC8irgWOBWgKq6EnjIbAcneSqwvapm/61YhCTrkmxIsgE4aJDXliRJQ/MqFlCvkCRJ02PJLTCAH1bVbTv16rh7juMfCzw9yYnAPsC+wOuB/ZLs0VphrAC2teO30ZtCbWuSPYAHAN/c+aJVdTa96dZoSQxJkjR+FlqvkCRJU2IQCYzPJfl1YPckq4DfA/7/2Q6uqlcArwBIchzwR1X1G0neA5wCXEBv5PGL2ynr2/Yn2/6PVFUNIG5JktQ9C6pXSJKGa67uFWAXCy2vQXQh+V3gYcCdwDuB2+jNMLJQLwdemmQzvTEuzmnl5wAHtvKXAmcsOWJJktRVg6pXSJKkCTOIFhi/UlVnAmfuKEjyDOA9uzqxqj4GfKytb6HX53XnY+4AnjGAOBfEzKIkSSOx6HqFJEmabINIYLyCn6xUzFTWGSYnJEnqrLGrV0iSfpIzeWg5LDqBkeQE4ETg0CRv6Nu1L3DXzGdJkiT9JOsVkiRpV5bSAuNrwAbg6UB/eu07wB8sJShJkjR1rFdIkqQ5LTqBUVVXAVcleWdV/XCAMUmSpCljvUKSJO3KIMbAOCzJXwFHAvvsKKyqhwzg2pIkabpYr5AkSTMaRALjH4BXAq8DHg+cxmCmZ5UkSdPHeoUkaVk4mcP4G0QC46eq6rIkqaqvAq9KshH47wO4tiRJmi7WKyRJGhPDTgoNIoFxZ5LdgOuTvBjYBtxvANeVJEnTx3qFpKFzyk9pPAyiSebvA/cBfg9YDfwmsHYA15UkSdPHeoUkSZrRkltgVNVn2urt9PqpSjMysy1J2hXrFZIkaTZLboGR5NIk+/Vt75/kQ0u9riRJmj7WKyRJ0mwG0YXkoKq6dcdGVd0CPHAA15UkSdPHeoUkSZrRIBIYdyf5mR0bSR4M1ACuK0mSpo/1CkmSNKNBzEJyJvCJJP8GBPgvwLoBXFeSJE0f6xWSJGlGgxjE81+SHAM8phW9pKpuXup1JUnS9LFeIUmSZrPoLiRJfr49HgP8DPC1tvxMK5MkSZoX6xWSJGlXltIC4w+B3wb+5wz7CnjCEq4tSZKmy6LrFUn2AT4O7E2vbnNRVb0yyeHABcCBwEbgOVX1gyR7A+cDq4FvAs+qqq8M8F4kSdIyWHQCo6p+uz0+fnDhSJKkabTEesWdwBOq6vYke9IbQ+ODwEuB11XVBUneApwOnNUeb6mqI5KcCrwWeNZAbkSSNBIbN26cdd/q1auHGImW06ITGEn+77n2V9X7FnttSZI0XZZSr6iqAm5vm3u2ZUerjV9v5ecBr6KXwDiprQNcBLwpSdp1JElSRy2lC8nT5thXgAkMSZI0X0uqVyTZnV43kSOANwNfAm6tqrvaIVuBQ9v6ocANAFV1V5Lb6HUzcbBQSZI6bCldSE4bZCCSJGl6LbVeUVU/Ao5Osh/wfuDnlxpTknX8eArXg5Z6PUmStDSLnoVkhyQHJnlDks8m2Zjk9UkOHERwkiRpuiy1XlFVtwIfBX4R2C/Jji9rVgDb2vo2YGX7eXsAD6A3mOfO1zq7qtZU1RpsnSFJ0sgtOYFBb3TvbwC/BpzS1t89gOtKkqTps+B6RZKfbi0vSPJTwJOA6+glMk5ph60FLm7r69s2bf9HHP9CkqTuW8oYGDscUlWv7tv+8ySO5D0G1qxZM+s+63GSpBFZTL3iEOC8Ng7GbsCFVfWBJNcCFyT5c+AK4Jx2/DnA25NsBr4FnDrYW5AkScthEAmMD7cpyC5s26cAHxrAdSVJ0vRZcL2iqjYBj5qhfAtw7AzldwDPWHqokiRpmAaRwPht4CXAP7bt3YDvJnkBvZnN9h3Az9CUcP5mSZp61iskSdKMlpzAqKr7DyIQSZIk6xWSJGk2S05gJPmlmcqr6uNLvbYkSZou1iskSTuzlbZ2GEQXkj/uW9+HXl/TjcATBnBtSZI0XaxXSJKkGQ2iC8nT+reTrAT+bqnXlSRJ08d6hSRJms1uy3DNrcAvzLYzyT5JPp3kqiSfS/KnrfzwJJ9KsjnJu5Ps1cr3btub2/7DliFmSZLUTXPWKyRJ0vQYxBgYbwSqbe4GHA18do5T7gSeUFW3J9kT+ESSDwIvBV5XVRckeQtwOnBWe7ylqo5o06q9FtjVfPBTbc2aNXPur6o590uSNCqLqFdIkqQpMYgxMDb0rd8FvKuq/n22g6v36fn2trlnW4pe39Zfb+XnAa+il8A4qa0DXAS8KUnKT+GSJE2iBdUrJEnS9BjEGBjnJfnptv6N+ZyTZHd6A3IdAbwZ+BJwa1Xd1Q7ZChza1g8FbmjXvyvJbcCBwM1LjV2SJHXLYuoVkqTxNNfsIuAMI/pJix4DIz2vSnIz8AXgi0m+keS/7+rcqvpRVR0NrKA3uvjPLzaOvnjWJdmQZANw0FKvJ0mShmcp9QpJkjQdljKI5x8AjwX+c1UdUFX7A48GHpvkD+Zzgaq6Ffgo8IvAfkl2tAhZAWxr69uAlQBt/wOAb85wrbOrak1VrcHWGZIkjZsl1yskSdJkW0oC4znAs6vqyzsKqmoL8JvAc2c7KclPJ9mvrf8U8CTgOnqJjFPaYWuBi9v6+rZN2/8Rx7+QJGniLKpeIUmSpsdSxsDYs6p+oqVDVX2jzS4ym0OA89o4GLsBF1bVB5JcC1yQ5M+BK4Bz2vHnAG9Pshn4FnDqEmKWJEndtNh6hSRJmhJLSWD8YDH7qmoT8KgZyrfQGw9j5/I7gGcsJkAt3VxTstoQRpI0QIuqV0iSpOmxlATGI5N8e4byAPss4bqSJGn6WK+QJElzWnQCo6p2H2QgkiRpelmvkCRJu7KUFhiSJEmSpGVid27p3pYyC4kkSZIkSdJQ2AJDY2fjxo1z7l+9evWQIpEkSZIk7cpcn+EW8vnNFhiSJEmSJKnzTGBIkiRJkqTOM4EhSZIkSZI6zwSGJEmSJEnqPAfxlCRJkjSRujr4+1zTo4JTpEqzsQWGJEmSJEnqPBMYkiRJkiSp80xgSJIkSZKkzjOBIUmSJEmSOs8EhiRJkiRJ6jwTGJIkSZIkqfNMYEiSJEmSpM7bY9QBSJIkSZIWZ82aNXPur6ohRSItP1tgSJIkSZKkzjOBIUmSJEmSOs8uJJpYGzdunHXf6tWrhxiJJEmSJGmpbIEhSZLGWpKVST6a5Nokn0vy+638gCSXJrm+Pe7fypPkDUk2J9mU5JjR3oEkSZoPExiSJGnc3QX8YVUdCTwGeFGSI4EzgMuqahVwWdsGOAFY1ZZ1wFnDD1mSJC2UXUgkSdJYq6obgRvb+neSXAccCpwEHNcOOw/4GPDyVn5+9YbmvzzJfkkOadeRNCbm6i4MdhmWJpEtMCRJ0sRIchjwKOBTwMF9SYmvAwe39UOBG/pO29rKJElSh9kCQ5IkTYQk9wPeC7ykqr6d5J59VVVJaoHXW0eviwnAQQMLVJIkLYoJDEmSNPaS7EkvefGOqnpfK75pR9eQJIcA21v5NmBl3+krWtm9VNXZwNnt+huWLXhJWmZr1qyZc3+vR53UfXYhkSRJYy29phbnANdV1d/27VoPrG3ra4GL+8qf22YjeQxwm+NfSJLUfbbAkDQvcw2U5SBZkkbsscBzgKuTXNnK/gR4DXBhktOBrwLPbPsuAU4ENgPfA04bbriSJGkxTGBIkqSxVlWfADLL7uNnOL6AFy1rUJIkaeCGnsBIshI4n95I4AWcXVWvT3IA8G7gMOArwDOr6pbWLPT19L4p+R7wvKr67LDj1uzsU6flMOktPib9/iRJmlZz1Y2tF0tLM4oWGHcBf1hVn01yf2BjkkuB5wGXVdVrkpwBnEFvrvYTgFVteTRwVnuUNADOoS5JkobFeoekpRh6AqMNknVjW/9Okuvozb1+EnBcO+w84GP0EhgnAee35p6XJ9lvx4jiw469C8zoSpIkSZKm0UjHwEhyGPAo4FPAwX1Jia/T62ICveTGDX2nbW1lU5nA0GANuxm/3QYGy29xevy9kiRJ0jQYWQIjyf3ozdf+kqr6dm+oi56qqiQLak6QZB2wrm0eNLBAJY01kxySJEnSZBhJAiPJnvSSF++oqve14pt2dA1JcgiwvZVvA1b2nb6ild1LVZ0NnN2uv2HZgpekITH5IkmSJP3YKGYhCXAOcF1V/W3frvXAWnpztq8FLu4rf3GSC+gN3nnbtI5/IU0CP5RLkiRJWoxRtMB4LPAc4OokV7ayP6GXuLgwyenAV4Fntn2X0JtCdTO9aVRPG2640nD5AX80BjmOhGNSSJIkSYM3illIPgFklt3Hz3B8AS9a1qAkSZIkSVKnjXQWEmkczLdFhDOaTD5bx0iSJEmjYwJD0sD4AV+SJA2CX9RImsluow5AkiRJkiRpV2yBIUkaOL85kyRpMq1Zs2bO/b0hDKXlYQsMSZIkSZLUebbAkCQtiK0rJEmSNAomMCRJkiTNykG6JXWFXUgkSZIkSVLn2QJDkqaA3T4kSZI07kxgSBPMJp+SJGlYTJZLWm4mMCRJkiRpF5w+VBo9ExiSJMAWO5IkSeo2ExiSJEnShDEprcWYq5WJLUzUBc5CIkmSJEmSOs8EhiRJkiRJ6jy7kKhTbLYmSZIkSZqJCQwNjckJSZIkSdJi2YVEkiRJkiR1ngkMSZIkSZLUeXYhkSSNhFP8SZIkaSFsgSFJkiRJkjrPBIYkSZIkSeo8u5BIkiRJU2qu7nx25dNSOAOhloMJDEmSJGmMOIaQpGllFxJJkiRJktR5tsCQJEljLcm5wFOB7VX18FZ2APBu4DDgK8Azq+qWJAFeD5wIfA94XlV9dhRxS+oOuztI48EWGJIkady9DXjKTmVnAJdV1SrgsrYNcAKwqi3rgLOGFKMkSVoiExiSJGmsVdXHgW/tVHwScF5bPw84ua/8/Oq5HNgvySHDiVSSJC2FXUgkSWPPUfQ1g4Or6sa2/nXg4LZ+KHBD33FbW9mNSJKkTjOBIUmSJlpVVZIFd2JPso5eNxOAgwYblSRJWqiRdCFJcm6S7Umu6Ss7IMmlSa5vj/u38iR5Q5LNSTYlOWYUMUuSpLFy046uIe1xeyvfBqzsO25FK/sJVXX5lZmxAAAgAElEQVR2Va2pqjXAzcsZrCRJ2rVRjYHxNhxsS5IkLZ/1wNq2vha4uK/8ue0LkscAt/V1NZEkSR02kgSGg21JkqRBSfIu4JPAQ5NsTXI68BrgSUmuB57YtgEuAbYAm4G/B35nBCFLkqRF6NIYGA62JUmSFqyqnj3LruNnOLaAFy1vRJIkaTl0KYFxj8UMtuVAW9NjzZo1c+7v1U0l6d7mmqkEnK1E0vLyPUiSlm5UY2DMZEmDbTnQliRJkiRJk6tLCQwH25IkSZIkSTMaSReSNtjWccBBSbYCr6Q3uNaFbeCtrwLPbIdfApxIb7Ct7wGnDT1gSdLIzNXs2ibXkiSNL7uGa6FGksBwsC1JkiRJkrQQnRzEU5KkLrNViKTl4vuLJM2uS2NgSJIkSZIkzcgEhiRJkiRJ6jy7kEiS1Mfm25IWyvcNafo4AOlomMCQJGkZzPWBBvxQI0mStFB2IZEkSZIkSZ1nC4wOmasZkk2QFs7nU5IkLYUtqSSpW0xgSJIkSZI6a9hfTPpFaHeZwNBU881J0qjNZ/A/vwWWJGnpHHhz/DkGhiRJkiRJ6jxbYEiSJGksDLI1klOfSpOlqy2rBxWXrUd6TGBIkjRF/NAmSZLGlQkMSZImgONkSJKkSWcCQ5IkSRPDZJ762exemiwmMCRJkjRydm+SJO2Ks5BIkiRJkqTOswWGJEm6F5vgS5KkLjKBIUmSNAFMPEnSeBrkWC1dnU52UExgSJIkadmYWJEkDYoJDGlAhp3tnPTsqiRpdEw6SJK6yASGtAuTPv3WpN+fJEmSpB8bVpeVhV5rPkxgSEM0n1YTXU0o2OJDkpbHsKcPHeTPc+pTSdIwmcCQJlhXkyGSJEmStFAmMCQNjAkTSVoejkkh/STrHdL02W3UAUiSJEmSJO2KLTAkdZJjbkgad7aakCSNs0HWxwd1LRMYkoaqq809u/gGLUmSJOnHTGBIGlvzSYaM8zRR82XCRJIkSdPABMYQdPVDjyRJkiRJ48IEhiSNwKS3mpj0+5PmGt/CsS0kSVoeY5PASPIU4PXA7sBbq+o1Iw5JkpbVsLq/TMsYH8NuDTffnzeuz+e4s14hjYbvjZKWYiwSGEl2B94MPAnYCnwmyfqquna0kfnmKmly+AGfef883/vHW5frFZIkaXZjkcAAjgU2V9UWgCQXACcBi6poDHvgP0mSlmKQ/7eWmqDpP26MDbReIXXRFPwdS5pCGYc3rySnAE+pqt9q288BHl1VL+47Zh2wrm3uU1UPH36kkiSp6+ZTr2jl1i0kSeqQcWmBsUtVdTZw9qjjkCRJk8G6hSRJ3TIuCYxtwMq+7RWtbE5JrgHuWK6gNKeDgJtHHcQU8nkfHZ/70fB5H41xb42wqHrFYkxZXWSa/h6918nkvU6mabrXcf//vEvjksD4DLAqyeH0KhinAr8+j/PuqKq5OwBqWSTZ4HM/fD7vo+NzPxo+76ORZMOoY1iixdYrFmNq6iLT9PfovU4m73UyTdu9jjqG5TYWCYyquivJi4EP0Zvu7Nyq+tyIw5IkSWPIeoUkSeNpLBIYAFV1CXDJqOOQJEnjz3qFJEnjZ7dRB7DMHHhrdHzuR8PnfXR87kfD5300fN7nb5qeK+91Mnmvk8l7nUwTf69jMY2qJEmSJEmabpPeAkOSJEmSJE0AExiSJEmSJKnzJjaBkeQpSb6QZHOSM0Ydz7RIcm6S7W3eew1JkpVJPprk2iSfS/L7o45pGiTZJ8mnk1zVnvc/HXVM0yTJ7kmuSPKBUccyTZJ8JcnVSa6chuna5iPJM9p7wN1JZpyqb6736SSvSrKtPadXJjlxeNEvzHzutR03Yz0syeFJPtXK351kr+FEvnBJDkhyaZLr2+P+Mxzz+L7X7cokdyQ5ue17W5Iv9+07evh3MT/zudd23I/67md9X/mkva5HJ/lk+13flORZffs6/7ru6nNQkr3b67S5vW6H9e17RSv/QpJfHmbcizGPe31pe9/dlOSyJA/u2zfj73NXzeNen5fkG3339Ft9+9a23/nrk6wdbuQDVlUTt9CbEu1LwEOAvYCrgCNHHdc0LMAvAccA14w6lmlagEOAY9r6/YEv+js/lOc9wP3a+p7Ap4DHjDquaVmAlwLvBD4w6limaQG+Ahw06ji6tAC/ADwU+BiwZpZjZn2fBl4F/NGo72OA9zprPQy4EDi1rb8F+G+jvqc57vV/AGe09TOA1+7i+AOAbwH3adtvA04Z9X0M8l6B22cpn6jXFfg5YFVb/0/AjcB+4/C6zudzEPA7wFva+qnAu9v6ke34vYHD23V2H/U9LfFeH9/3N/nfdtxr257x97mLyzzv9XnAm2Y49wBgS3vcv63vP+p7WuwyqS0wjgU2V9WWqvoBcAFw0ohjmgpV9XF6/7w1RFV1Y1V9tq1/B7gOOHS0UU2+6rm9be7ZFkdGHoIkK4BfAd466likqrquqr6wi2Mm4n16PvfKLPWwJAGeAFzUjjsPOHn5ol2yk+jFCPOL9RTgg1X1vWWNanks9F7vMYmva1V9saqub+tfA7YDPz20CJdmPp+D+p+Di4Dj2+t4EnBBVd1ZVV8GNrfrddUu77WqPtr3N3k5sGLIMQ7KUj7f/jJwaVV9q6puAS4FnrJMcS67SU1gHArc0Le9lTGsJEiL0ZoBPopeawAts/S6MVxJr3JzaVX5vA/H3wEvA+4edSBTqIAPJ9mYZN2ogxlHs7xPv7g1cT53tub7Y2S2etiBwK1VdddO5V11cFXd2Na/Dhy8i+NPBd61U9lftNf1dUn2HniEgzPfe90nyYYkl+/oKsOEv65JjqX3jfeX+oq7/LrO53PQPce01+02eq/juH2GWmi8pwMf7Nue6fe5q+Z7r7/WfjcvSrJygeeOhT1GHYCkwUlyP+C9wEuq6tujjmcaVNWPgKOT7Ae8P8nDq8oxYJZRkqcC26tqY5LjRh3PFHpcVW1L8kDg0iSfb63vJlqSfwUeNMOuM6vq4gVcZ6b36bOAV9NLDr0a+J/A85cW8eIN6l7HwVz32r9RVZVk1hZ2SQ4BjgI+1Ff8CnofkPcCzgZeDvzZUmNerAHd64Pb3/9DgI8kuZreh99OGfDr+nZgbVXtSJh36nXV/CT5TWAN8F/7in/i97mqvjTzFcbCPwHvqqo7k7yAXiubJ4w4poGb1ATGNmBl3/aKViZNrCR70qsUv6Oq3jfqeKZNVd2a5KP0muSZwFhejwWent5Ah/sA+yb5x6r6zRHHNRWqalt73J7k/fSatU58AqOqnrjUa8z2Pl1VN/Ud8/fASAemHcC9zlYP+yawX5I92re+I6+fzXWvSW5KckhV3dg+yG6f41LPBN5fVT/su/aOb/nvTPIPwB8NJOhFGsS99v39b0nyMXotid7LBL6uSfYF/ple4u7yvmt36nWdwXw+B+04ZmuSPYAH0Pv7HLfPUPOKN8kT6SWv/mtV3bmjfJbf564mMHZ5r1X1zb7Nt9Ib72XHucftdO7HBh7hkExqF5LPAKvSGxF5L3pN+jo/sqy0WK3f4jnAdVX1t6OOZ1ok+enW8oIkPwU8Cfj8aKOafFX1iqpaUVWH0Xt//4jJi+FIct8k99+xDjwZE3bzMtf7dPsQtcOvMv7P6Yz1sKoq4KP0xooAWAt0uUXHenoxwq5jfTY7dR/Z8bq21/5kuv267vJek+y/o7tEkoPoJZOvncTXtf3evh84v6ou2mlf11/X+XwO6n8OTqH3f7Ra+anpzVJyOLAK+PSQ4l6MXd5rkkcB/wt4elVt7yuf8fd5aJEv3Hzutf9/ydPpjbUEvZZhT273vD+9/939rcXGyyBHBO3SApxIb4TvL9HLnI48pmlY6P3zvhH4Ib3+VaePOqZpWIDH0Wt6vAm4si0njjquSV+ARwBXtOf9GuC/jzqmaVvofaPgLCTDe74fQm/k86uAz/n/9Z7n5Vfb/7w7gZuAD7Xy/wRc0tZnfZ+m10T96rZvPXDIqO9pKffatmesh7XfoU/TGxzwPcDeo76nOe71QOAy4HrgX4EDWvka4K19xx1G7xvO3XY6/yPtdb0G+EfarFVdXOZzr8D/1e7nqvZ4et/5E/W6Ar/Z6rJX9i1Hj8vrOtPfH71uLk9v6/u012lze90e0nfume28LwAnjPpeBnCv/9req3a8jut39fvc1WUe9/pX9P43X0Uvqfjzfec+v73em4HTRn0vS1nSbkiSJEmSJKmzJrULiSRJkiRJmiAmMCRJkiRJUueZwJAkSZIkSZ1nAkOSJEmSJHWeCQxJkiRJktR5JjAkSZIkSVLnmcCQJEmSJEmdZwJD0kgkeUGSt7T1PZO8Pcl5SfYcdWySJGl4kty+hHPPTbI9yTU7lT8lyReSbE5yRl/5TyX5tyS7t+1510eS7JXk40n2WGy8kpbGBIakUTkK2JRkX+CDwH9U1dqq+uGI45IkSePjbcBT+gtacuLNwAnAkcCzkxzZdj8feF9V/ahtz7s+UlU/AC4DnrUcNyJp10xgSBqVRwDfAj4GvKeqzhxtOJIkaZSSvDTJNW15SV/5/9NaU3wiybuS/NGOfVX1cXr1iX7HApuraktLOlwAnNT2/QZwcd+xC62P/O92DUkjYPMnSaNyFPBG4PlV9U+jDkaSJI1OktXAacCjgQCfSvJv9D6v/BrwSGBP4LPAxl1c7lDghr7trcCjk+wFPKSqvtK3b6H1kWuA/zyP4yQtAxMYkoYuyUrgduB64JC+8ocBZwI3A1ur6n+MJkJJkjRkjwPeX1XfBUjyPuC/0GsxfnFV3QHckWQpX3ocBNy6Y2O2+kjbN2OdpKp+lOQHSe5fVd9ZQiySFsEEhqRROAq4Cvht4PIkn6mqK4AnA2+vqg+ONDpJkjTOtgEr+7ZXtLLvA/v0lc9WH4G56yR7A3cMPGpJu+QYGJJG4RHA1VV1I/BbwLuTPAA4B3hcknOSvGCkEUqSpGH6/4CTk9wnyX2BX21l/w48Lck+Se4HPHUe1/oMsCrJ4a3byKnA+qq6Bdg9yY4kxmz1EZilTpLkQOBmBx2XRsMEhqRROAq4GqCqLgUuBM6tqm9X1ZlVdTrw9CS+R0mSNAWq6rP0ZhT5NPAp4K1VdUVVfQZYD2yiN0vI1cBtO85L8i7gk8BDk2xNcnpV3QW8GPgQcB1wYVV9rp3yYXrdVWCW+kjbnq1O8njgn5fhKZA0D6mqUccgSQAkORn4ZeAu4I6q+uMRhyRJkkYsyf2q6vYk9wE+DqxrCY/FXOsY4A+q6jm7OG7GOkkbm+OMqvriYn6+pKUxgSFJkiSps5K8EziS3vgV51XVXy3xes9v1/nRAs/bCzi1qs5fys+XtHgmMCRJkiRJUufZv1ySJEmSJHWeCQxJkiRJktR5JjAkSZIkSVLnmcCQJEmSJEmdZwJDkiRJkiR1ngkMSZIkSZLUeSYwJEmSJElS55nAkCRJkiRJnWcCQ5IkSZIkdZ4JDEmSJEmS1HkmMCRJkiRJUueZwJAkSZIkSZ1nAkOSJEmSJHWeCQxJkiRJktR5JjAkSZIkSVLn7THqAJZZjToASZImUEYdwAhZt5AkafDmVbewBYYkSZIkSeo8ExiSJEmSJKnzTGBIkiRJkqTOM4EhSZIkSZI6zwSGJEmSJEnqPBMYkiRJkiSp80xgSJIkSZKkzjOBIUmSJEmSOs8EhiRJkiRJ6jwTGJIkSZIkqfNMYEiSpLGQ5Nwk25NcM8O+P0xSSQ5q20nyhiSbk2xKckzfsWuTXN+WtcO8B0mStHgmMCRJ0rh4G/CUnQuTrASeDPxHX/EJwKq2rAPOasceALwSeDRwLPDKJPsva9SSJGkg9hh1AJI0jjZu3DjrvtWrVw8xEml6VNXHkxw2w67XAS8DLu4rOwk4v6oKuDzJfkkOAY4DLq2qbwEkuZReUuRdyxi6pEXwf62kndkCQ5Ikja0kJwHbquqqnXYdCtzQt721lc1WLkmSOs4WGJIkaSwluQ/wJ/S6jyzH9dfR634CsE9VPXw5fo4kSZqfZUtgJDkXeCqwfcc//CR/DTwN+AHwJeC0qrq17XsFcDrwI+D3qupDrfwpwOuB3YG3VtVrlitmSRokm75Ky+5ngcOBq5IArAA+m+RYYBuwsu/YFa1sG71uJP3lH5vp4lV1NnA2QJINgw1dkiQt1HJ2IXkbPznQ1qXAw6vqEcAXgVcAJDkSOBV4WDvn/02ye5LdgTfTG4jrSODZ7VhJkjTlqurqqnpgVR1WVYfR6w5yTFV9HVgPPLfNRvIY4LaquhH4EPDkJPu3wTuf3MokSVLHLVsCo6o+Dnxrp7IPV9VdbfNyet96QG+grQuq6s6q+jKwmd7I4McCm6tqS1X9ALigHStJkqZMkncBnwQemmRrktPnOPwSYAu9OsXfA78D0AbvfDXwmbb82Y4BPSVJUreNcgyM5wPvbuuH0kto7NA/oNbOA209evlDkyRJXVNVz97F/sP61gt40SzHnQucO9DgJEnSshvJLCRJzgTuAt4xwGuuS7KhLdcM6rqSJEmSJGn0ht4CI8nz6A3ueXz7dgRmH2iLOcrvxYG2JEmSpPHgQNeSFmOoLTDajCIvA55eVd/r27UeODXJ3kkOB1YBn6bXN3VVksOT7EVvoM/1w4xZkiRJkiSN3nJOo/ouetOUHZRkK/BKerOO7A1c2qY7u7yqXlhVn0tyIXAtva4lL6qqH7XrvJje6OC7A+dW1eeWK2ZJ2sFvhiRJkqRuWbYExiwDbZ0zx/F/AfzFDOWX0BtJXJIkSZIkTamRDOIpSZIkSZK0ECYwJEmSJElS55nAkCRJkiRJnTf0aVQlST/mYKGSJEnS/NgCQ5IkSZIkdZ4JDEmSJEmS1HkmMCRJkiRJUueZwJAkSZIkSZ1nAkOSJEmSJHWeCQxJkiRJktR5JjAkSZIkSVLnmcCQJEmSJEmdZwJDkiRJkiR13h6jDkCSJEnS5Ni4ceOs+1avXj3ESCRNGltgSJIkSZKkzjOBIUmSJEmSOs8EhiRJkiRJ6jzHwJCkDrMfsSRJktRjCwxJkiRJktR5JjAkSZIkSVLnmcCQJEmSJEmdZwJDkiRJkiR1ngkMSZIkSZLUeSYwJEmSJElS55nAkCRJYyHJuUm2J7mmr+yvk3w+yaYk70+yX9++VyTZnP/T3r1HWVqVdx7//qS5iDeQnhDS3QYcO2aIV6oDZDHJUtoo4KWJQQZipEWSzmTw7kQxzARXjGt0JROCxjGrR4htFnIJktBjiMighEkmELoQkUuUHozSHRCRixpitOMzf5zdSVl2dVdVn8t7Tn0/a9Wq9917n/c89XZVvbuf2pfkC0leMqP8hFa2Nck5w/46JEnS4pjAkCRJ4+IjwAmzyq4FnlVVzwG+CLwTIMmRwGnAT7TX/I8k+yTZB/ggcCJwJHB6aytJkjpu2agDkCRJmo+quiHJ4bPKPjXj9EbglHa8Dri0qv4J+FKSrcDRrW5rVd0DkOTS1vbOAYYuaQCmp6fnrJuamhpiJJKGxREYkiRpUrwO+PN2vAK4d0bdtlY2V/kPSLIhyZYkW4Dl/Q9XkiQthAkMSZI09pKcC+wALu7XNatqY1Wtqao1wIP9uq4kSVocp5BIWnIccipNliSvBV4GrK2qasXbgVUzmq1sZeymXJIkdZgjMCRJ0thKcgLwduAVVfXYjKrNwGlJ9k9yBLAa+BvgZmB1kiOS7Edvoc/Nw45bkiQtnCMwJEnSWEhyCfACYHmSbcB59HYd2R+4NgnAjVX1H6vqjiSX01uccwdwdlX9c7vO64FrgH2Ai6rqjqF/MZIkacFMYEiSpLFQVafvovjC3bR/D/CeXZRfDVzdx9AkSdIQOIVEkiRJkiR13sASGEkuSvJAkttnlD01ybVJ7m6fD27lSfL+JFuT3JbkqBmvWd/a351k/aDilSRJkiRJ3TXIERgfAU6YVXYOcF1VrQaua+cAJ9JbXGs1sAH4EPQSHvTmtx4DHA2ctzPpIUmSJEmSlo6BJTCq6gbgoVnF64BN7XgTcPKM8o9Wz43AQUkOA14CXFtVD1XVw8C1/GBSRJIkSZIkTbhhr4FxaFXd147vBw5txyuAe2e029bK5iqXJEmSJElLyMh2IamqSlL9ul6SDfSmnwAs79d1JUmSJEnS6A17BMZX29QQ2ucHWvl2YNWMditb2VzlP6CqNlbVmqpaAzzY78AlSZIkSdLoDDuBsRnYuZPIeuCqGeVntN1IjgUebVNNrgFenOTgtnjni1uZJEmSJElaQgY2hSTJJcALgOVJttHbTeS9wOVJzgK+DJzaml8NnARsBR4DzgSoqoeSvBu4ubX7zaqavTCoJEmSJEmacANLYFTV6XNUrd1F2wLOnuM6FwEX9TE0SZIkSZI0ZoY9hUSSJEmSJGnBTGBIkiRJkqTOM4EhSZIkSZI6zwSGJEmSJEnqPBMYkiRJkiSp80xgSJIkSZKkzjOBIUmSJEmSOs8EhiRJkiRJ6jwTGJIkSZIkqfNMYEiSJEmSpM5bNuoAJEmSJI2P6enpOeumpqaGGImkpcYRGJIkSZIkqfNMYEiSJEmSpM4zgSFJkiRJkjrPBIYkSZIkSeo8ExiSJEmSJKnzTGBIkiRJkqTOM4EhSZIkSZI6zwSGJEkaC0kuSvJAkttnlD01ybVJ7m6fD27lSfL+JFuT3JbkqBmvWd/a351k/Si+FkmStHAmMCRJ0rj4CHDCrLJzgOuqajVwXTsHOBFY3T42AB+CXsIDOA84BjgaOG9n0kOSJHWbCQxJkjQWquoG4KFZxeuATe14E3DyjPKPVs+NwEFJDgNeAlxbVQ9V1cPAtfxgUkSSJHWQCQxJkjTODq2q+9rx/cCh7XgFcO+Mdtta2VzlkiSp45aNOgBJkqR+qKpKUv26XpIN9KafACzv13UlSdLimMCQJEnj7KtJDquq+9oUkQda+XZg1Yx2K1vZduAFs8qv39WFq2ojsBEgyZb+hi1pGKanp+esm5qaGmIkkvrBKSSSJGmcbQZ27iSyHrhqRvkZbTeSY4FH21STa4AXJzm4Ld754lYmSZI6zhEYkiRpLCS5hN7oieVJttHbTeS9wOVJzgK+DJzaml8NnARsBR4DzgSoqoeSvBu4ubX7zaqavTCoJEnqIBMYkiRpLFTV6XNUrd1F2wLOnuM6FwEX9TE0SZI0BE4hkSRJkiRJnWcCQ5IkSZIkdZ4JDEmSJEmS1HkmMCRJkiRJUueZwJAkSZIkSZ3nLiSSNOamp6fnrJuamhpiJJIkSdLgmMCQNHHm+g+9/5mXJEmSxtdIppAkeUuSO5LcnuSSJAckOSLJTUm2JrksyX6t7f7tfGurP3wUMUuSJEmSpNEZegIjyQrgjcCaqnoWsA9wGvA+4PyqegbwMHBWe8lZwMOt/PzWTpIkSZIkLSGjWsRzGfD4JMuAA4H7gOOBK1r9JuDkdryundPq1ybJEGOVJEmSJEkjNvQERlVtB34H+Aq9xMWjwDTwSFXtaM22ASva8Qrg3vbaHa39IbOvm2RDki1JtgDLB/pFSJIkSZKkoRrFFJKD6Y2qOAL4EeAJwAl7e92q2lhVa6pqDfDg3l5PkiRJkiR1xyimkLwI+FJVfa2qvgtcCRwHHNSmlACsBLa34+3AKoBW/xTg68MNWZIkSZIkjdIoEhhfAY5NcmBby2ItcCfwGeCU1mY9cFU73tzOafWfrqoaYrySJEmSJGnERrEGxk30FuO8Bfh8i2Ej8A7grUm20lvj4sL2kguBQ1r5W4Fzhh2zJEmSJEkarWV7btJ/VXUecN6s4nuAo3fR9tvAq4YRlyRJkiRJ6qZRbaMqSZIkSZI0byYwJEmSJElS55nAkCRJkiRJnTevBEaSf5tk/3b8giRvTHLQYEOTJEmTyr6FJElaqPmOwPg48M9JnkFvx5BVwMcGFpUkSZp09i0kSdKCzDeB8b2q2gH8HPCBqvo14LDBhSVJkiacfQtJkrQg801gfDfJ6cB64BOtbN/BhCRJkpYA+xaSJGlB5pvAOBP4KeA9VfWlJEcAfzS4sCRJ0oSzbyFJkhZk2XwaVdWdSd4BPK2dfwl43yADkyRJk8u+hSRJWqj57kLycuBW4JPt/HlJNg8yMEmSNLnsW0iSpIWa1wgM4F3A0cD1AFV1a5KnDygmSZI0+d6FfQupc6anp+esm5qaGmIkkvSD5r2IZ1U9Oqvse/0ORpIkLRn2LSRJ0oLMN4FxR5JfAPZJsjrJB4D/O8C4JEnSZOtr3yLJW5LckeT2JJckOSDJEUluSrI1yWVJ9mtt92/nW1v94f35kiRJ0iDNN4HxBuAngH8CPgY8CrxpUEFJkqSJ17e+RZIVwBuBNVX1LGAf4DR6i4KeX1XPAB4GzmovOQt4uJWfj4uHSpI0FuabwHhpVZ1bVT/ZPv4L8IpBBiZJkiZav/sWy4DHJ1kGHAjcBxwPXNHqNwEnt+N17ZxWvzZJ9uK9JUnSEMw3gfHOeZZJkiTNR9/6FlW1Hfgd4Cv0EhePAtPAI1W1ozXbBqxoxyuAe9trd7T2hyzmvSVJ0vDsdheSJCcCJwErkrx/RtWTgR27fpUkSdKuDaJvkeRgeqMqjgAeAf4YOGEvQyXJBmBDO12+t9eTJEl7Z0/bqP49sIXekM6Zeyp9E3jLoIKSJEkTaxB9ixcBX6qqrwEkuRI4DjgoybI2ymIlsL213w6sAra1KSdPAb4++6JVtRHY2K65ZZGxSZKkPtltAqOqPgd8LsnHquq7Q4pJkiRNqAH1Lb4CHJvkQOAfgbX0kiSfAU4BLgXWA1e19pvb+V+3+k9XVfUpFkmSNCB7GoGx0+FJ/htwJHDAzsKqevpAopIkSZOub32LqropyRXALfSmoXyW3siJPwMuTfJbrezC9pILgT9KshV4iN6OJZIkqePmm8D4Q+A8eluNvRA4k/kvACpJkjRbX/sWVXVeu95M9wBH76Ltt4FXLfa9JEnSaMy3o/D4qrHmgAwAABifSURBVLoOSFV9uareBbx0cGFJkqQJZ99CkiQtyHxHYPxTkscBdyd5Pb3Fr544uLAkSdKEs28hSZIWZL4jMN4EHAi8EZgCfpHe4leSJEmLYd9CkiQtyLxGYFTVze3wW/TmqEqSJC2afQtJkrRQ80pgJLkWeFVVPdLODwYuraqXDDI4SZptenp6zrqpqakhRiJpb9i3kCRJCzXfKSTLd3YwAKrqYeCHBhOSJElaAuxbSJKkBZlvAuN7SZ628yTJjwI1mJAkSdISYN9CkiQtyHx3ITkX+MskfwEE+Glgw8CikiRJk86+haROmGt6qlNTpe6Z7yKen0xyFHBsK3pzVT04uLAkSdIks28hSZIWardTSJL8ePt8FPA04O/bx9NamSRJ0rzZt5AkSYu1pxEYbwN+Gfjvu6gr4Pi+RyRJkiaZfQtJkrQou01gVNUvt88vHE44kiRpktm3kCRJi7XbBEaSV+6uvqqu7G84kiRpktm3kCRJi7WnKSQv301dAYvqZCQ5CPgw8Kx2ndcBXwAuAw4H/g44taoeThLgAuAk4DHgtVV1y2LeV5IkjdxA+haSJGny7WkKyZkDet8LgE9W1SlJ9gMOBH4duK6q3pvkHOAc4B3AicDq9nEM8KH2WZIkjZkB9i0kSdKE2+0uJDslOSTJ+5PckmQ6yQVJDlnMGyZ5CvAzwIUAVfWdqnoEWAdsas02ASe343XAR6vnRuCgJIct5r0lSVI39LNvIUmSloZ5JTCAS4GvAT8PnNKOL1vkex7RXv+HST6b5MNJngAcWlX3tTb3A4e24xXAvTNev62VfZ8kG5JsSbIFWL7I2CRJ0nD0s28hSZKWgPkmMA6rqndX1Zfax2/xrwmGhVoGHAV8qKqeD/wDveki/6Kqit482Hmrqo1Vtaaq1gAPLjI2SZI0HP3sW0iSpCVgvgmMTyU5Lcnj2sepwDWLfM9twLaquqmdX0EvofHVnVND2ucHWv12YNWM169sZZIkaXz1s28hSZKWgPkmMH4Z+BjwnfZxKfArSb6Z5BsLecOquh+4N8kzW9Fa4E5gM7C+la0HrmrHm4Ez0nMs8OiMqSaSJGk89a1vIUmSloY9baMKQFU9qc/v+wbg4rYDyT3AmfSSKZcnOQv4MnBqa3s1vS1Ut9LbRtXVyyVJGnMD6FtIkqQJN68ERpKf2VV5Vd2wmDetqluBNbuoWruLtgWcvZj3kST1TE9P77J8ampqyJFIPf3uW0iSpMk3rwQG8Gszjg8AjgamgeP7HpEkSVoK7FtIkqQFme8UkpfPPE+yCvi9gUQkSZImnn0LSZK0UPNdxHO2bcC/62cgkiRpSbNvIUmSdmu+a2B8AKh2+jjgecAtgwpKkiRNNvsWkiRpoea7BsaWGcc7gEuq6q8GEI8kSVoa7FtIkqQFme8aGJuS/Jt2/LXBhiRJkiadfQtJkrRQu10DIz3vSvIg8AXgi0m+luQ3hhOeJEmaJPYtJEnSYu1pEc+3AMcBP1lVT62qg4FjgOOSvGXg0UmSpEkzkL5FkoOSXJHkb5PcleSnkjw1ybVJ7m6fD25tk+T9SbYmuS3JUf350iRJ0iDtKYHxGuD0qvrSzoKqugf4ReCMQQYmSZIm0qD6FhcAn6yqHweeC9wFnANcV1WrgevaOcCJwOr2sQH40F68ryRJGpI9JTD2raoHZxe2uar7DiYkSZI0wfret0jyFOBngAvbtb5TVY8A64BNrdkm4OR2vA74aPXcCByU5LDFvLckSRqePSUwvrPIOkmSpF0ZRN/iCOBrwB8m+WySDyd5AnBoVd3X2twPHNqOVwD3znj9tlYmSZI6bE+7kDw3yTd2UR7ggAHEI0mSJtsg+hbLgKOAN1TVTUku4F+niwBQVZWkFnLRJBvoTTEBWL7I2CRJUp/sNoFRVfsMKxBJkjT5BtS32AZsq6qb2vkV9BIYX01yWFXd16aIPNDqtwOrZrx+ZSubHetGYCNAki0DiFuSJC3AnqaQSJIkdVpV3Q/cm+SZrWgtcCewGVjfytYDV7XjzcAZbTeSY4FHZ0w1kSRJHbWnKSSSJEnj4A3AxUn2A+4BzqT3h5rLk5wFfBk4tbW9GjgJ2Ao81tpKkqSOM4EhSZLGXlXdCqzZRdXaXbQt4OyBByVJkvrKKSSSJEmSJKnzTGBIkiRJkqTOM4EhSZIkSZI6zzUwJEmSpCVienp6zrqpqakhRiJJC+cIDEmSJEmS1HkmMCRJkiRJUueZwJAkSZIkSZ1nAkOSJEmSJHWeCQxJkiRJktR5JjAkSZIkSVLnmcCQJEmSJEmdZwJDkiRJkiR13rJRByBJO01PT89ZNzU1NcRIJEmSJHWNIzAkSZIkSVLnmcCQJEmSJEmd5xQSSZIkSZqDU1yl7nAEhiRJkiRJ6ryRJTCS7JPks0k+0c6PSHJTkq1JLkuyXyvfv51vbfWHjypmSZIkSZI0GqMcgfEm4K4Z5+8Dzq+qZwAPA2e18rOAh1v5+a2dJEmSJElaQkaSwEiyEngp8OF2HuB44IrWZBNwcjte185p9Wtbe0mSJEmStESMagTG7wFvB77Xzg8BHqmqHe18G7CiHa8A7gVo9Y+29pIkSZIkaYkYegIjycuAB6pq7uV8F3fdDUm2JNkCLO/ntSVJkiRJ0miNYhvV44BXJDkJOAB4MnABcFCSZW2UxUpge2u/HVgFbEuyDHgK8PXZF62qjcBGgJbEkCRJkiRJE2LoIzCq6p1VtbKqDgdOAz5dVa8GPgOc0pqtB65qx5vbOa3+01VVQwxZkiRJkiSN2Ch3IZntHcBbk2ylt8bFha38QuCQVv5W4JwRxSdJkiRJkkZkFFNI/kVVXQ9c347vAY7eRZtvA68aamCSJEmSJKlTujQCQ5IkSZIkaZdMYEiSJEmSpM4zgSFJkiZCkn2SfDbJJ9r5EUluSrI1yWVJ9mvl+7fzra3+8FHGLUmS5scEhiRJmhRvAu6acf4+4PyqegbwMHBWKz8LeLiVn9/aSZKkjjOBIUmSxl6SlcBLgQ+38wDHA1e0JpuAk9vxunZOq1/b2kuSpA4zgSFJkibB7wFvB77Xzg8BHqmqHe18G7CiHa8A7gVo9Y+29pIkqcNGuo2qJHXVmjVr5qyrqiFGImlPkrwMeKCqppO8oI/X3QBsaKfL+3VdSZK0OCYwJEnSuDsOeEWSk4ADgCcDFwAHJVnWRlmsBLa39tuBVcC2JMuApwBfn33RqtoIbARIsmXgX4UkSdotp5BIkqSxVlXvrKqVVXU4cBrw6ap6NfAZ4JTWbD1wVTve3M5p9Z8uh1ZJktR5JjAkSdKkegfw1iRb6a1xcWErvxA4pJW/FThnRPFJkqQFcAqJpCXH9S2kyVVV1wPXt+N7gKN30ebbwKuGGpg0BNPT03PWTU1NDTESSRoMR2BIkiRJkqTOcwSGJEmSpL6Zz0hHR0NKWgwTGJIkhx1LkiSp80xgSNIi+JcjSZIkabhcA0OSJEmSJHWeCQxJkiRJktR5JjAkSZIkSVLnuQaGpKFxoUhJkjSJ7ONIw+EIDEmSJEmS1HkmMCRJkiRJUueZwJAkSZIkSZ3nGhiSJEmSOmfNmjVz1lXVECOR1BWOwJAkSZIkSZ1nAkOSJEmSJHWeCQxJkiRJktR5JjAkSZIkSVLnuYinpIkz16JfLvglSdLec3FNSaPiCAxJkiRJktR5JjAkSZIkSVLnOYVE0lhx2KokSZK0NDkCQ5IkSZIkdZ4jMCRpQBwtIkmSJPWPIzAkSZIkSVLnDX0ERpJVwEeBQ4ECNlbVBUmeClwGHA78HXBqVT2cJMAFwEnAY8Brq+qWYcctafemp6fnrJuamhpiJJIkabEcPSipy0YxAmMH8LaqOhI4Fjg7yZHAOcB1VbUauK6dA5wIrG4fG4APDT9kSZIkSZI0SkNPYFTVfTtHUFTVN4G7gBXAOmBTa7YJOLkdrwM+Wj03AgclOWzIYUuSJEmSpBEa6RoYSQ4Hng/cBBxaVfe1qvvpTTGBXnLj3hkv29bKJEmSJEnSEjGyXUiSPBH4OPDmqvpGb6mLnqqqJAuaZJdkA70pJgDL+xaoJEmSNGKuNSVJI0pgJNmXXvLi4qq6shV/NclhVXVfmyLyQCvfDqya8fKVrez7VNVGYGO7/paBBS9JkiRp5FxwVFp6hj6FpO0qciFwV1X97oyqzcD6drweuGpG+RnpORZ4dMZUE0mStMQlWZXkM0nuTHJHkje18qcmuTbJ3e3zwa08Sd6fZGuS25IcNdqvQJIkzcco1sA4DngNcHySW9vHScB7gZ9NcjfwonYOcDVwD7AV+J/AfxpBzJIkqbvc4UySpCVg6FNIquovgcxRvXYX7Qs4e6BBSdKIOPxV2nttZOZ97fibSWbucPaC1mwTcD3wDmbscAbcmOSgndNYhx27JEmav5HuQiJJktRP/dzhLMmGJFva2louEC5J0oiZwJAkSRNh9g5nM+vaaIsFDWuqqo1Vtaaq1gAP9i9SSZK0GCYwJEnS2NvdDmetfsE7nEmSpG4xgSFJksaaO5xJkrQ0DH0RT0mSpD7bucPZ55Pc2sp+nd6OZpcnOQv4MnBqq7saOIneDmePAWcON1xJkrQYJjAkSdJYc4czSV03PT09Z93U1NQQI5HGmwkMSZIkaYlw+25J48wEhqQ9GtZfDexUSZIkSZqLi3hKkiRJkqTOM4EhSZIkSZI6zwSGJEmSJEnqPNfAkCRJkkbMXSokac8cgSFJkiRJkjrPERiS1GHuzCJJkiT1mMCQNDT+Z1ySJEnSYpnAkCRJkjSR/OOJNFlMYEiS5sUF5iSp2/zPuqRJ5yKekiRJkiSp8xyBIUmSJEkj5khHac9MYEjygSlJkiSp85xCIkmSJEmSOs8RGJIkSSM210g4R8FJkvSvTGBIkiRJA+RUTUnqD6eQSJIkSZKkznMEhjThhvVXH/eelyRpcHzOSpIjMCRJkiRJ0hgwgSFJkiRJkjrPKSSSJEmSlqy5puc4NUfqHhMY0phzZXN1hd+LkiQNls9aLXUmMCRJkqS9MNd/Kv0PpST1lwkMSXvkyueSpKXK3bwkqTtMYEgd5jBBzcd8Or3O7/1+/frZGtbPqH/dlcabyQlJ6g8TGJIkDYhJSKnb/BmVpPFiAkOSfxlSp+xptMGoRz30+30kDY4/x1pq/J4fHe/9cIxNAiPJCcAFwD7Ah6vqvSMOSQswbj/Qwxiu7ZxaLUUmH37QOMU6SexXSJqvSetL9eO5M59r7KlNP66xFC31vtRYJDCS7AN8EPhZYBtwc5LNVXXnKOLp6j/moEza19uVX9r9MmkPVUkatK71K7Q4PmclQXdGbs5HP5I6S91YJDCAo4GtVXUPQJJLgXWAHY0JMqzEwjix0yRpPibtd98QdK5f4UK7C9elZ6T/furS9+PeWorPlEn7mid58e9xSWCsAO6dcb4NOGZEsfRFl0YBDCOWfv1SH/VfW/a0o8PMNot9j4VcQ+qXvdnJZGcbv6c1RsauXzFOSf5+9D2SzFm/kN8n/t5SV/TjOavRWSoJDuj+15Nx+GFIcgpwQlX9Ujt/DXBMVb1+RpsNwIZ2ekBVPWv4kUqSpK6bT7+ildu3kCSpQ8ZlBMZ2YNWM85Wt7F9U1UZg4zCDkiRJY2mP/QqwbyFJUteMywiMZcAXgbX0Ohg3A79QVXeMIJbbgW8P+30nwHLgwVEHMaa8d4vjfVsc79vi+Nf5MdKlfkWLx77FYPj7bDC8r4PhfR0c7+1gjKTvMxYjMKpqR5LXA9fQ2+7solF1MoBvV9XcE9S0S0m2eN8Wx3u3ON63xfG+LU6SLaOOQfPXsX4F2LcYCH+fDYb3dTC8r4PjvR2MUfV9xiKBAVBVVwNXjzoOSZI0/uxXSJI0fh436gAkSZIkSZL2xATGwrmY1+J43xbPe7c43rfF8b4tjvdNe8Pvn8Hwvg6G93UwvK+D470djJHc17FYxFOSJEmSJC1tjsCQJEmSJEmdZwJjLyR5W5JKsnzUsYyDJO9OcluSW5N8KsmPjDqmcZDkt5P8bbt3f5LkoFHHNC6SvCrJHUm+l8TVp/cgyQlJvpBka5JzRh3POEhyUZIH2jaY0qL4fBwcn6GD4fO1v3z+DobP6P5LsirJZ5Lc2X4HvGnYMZjAWKQkq4AXA18ZdSxj5Ler6jlV9TzgE8BvjDqgMXEt8Kyqeg7wReCdI45nnNwOvBK4YdSBdF2SfYAPAicCRwKnJzlytFGNhY8AJ4w6CI09n4+D4zN0MHy+9onP34H6CD6j+20H8LaqOhI4Fjh72N+vJjAW73zg7YCLiMxTVX1jxukT8N7NS1V9qqp2tNMbgZWjjGecVNVdVfWFUccxJo4GtlbVPVX1HeBSYN2IY+q8qroBeGjUcWi8+XwcHJ+hg+Hzta98/g6Iz+j+q6r7quqWdvxN4C5gxTBjWDbMN5sUSdYB26vqc0lGHc5YSfIe4AzgUeCFIw5nHL0OuGzUQWgirQDunXG+DThmRLFIS47Px6HwGaou8vmrsZTkcOD5wE3DfF8TGHNI8r+BH95F1bnAr9ObPqJZdnffquqqqjoXODfJO4HXA+cNNcCO2tN9a23OpTds6+JhxtZ187l3kjRqPh8Hx2foYPh8lTSXJE8EPg68edYowoEzgTGHqnrRrsqTPBs4Atg5+mIlcEuSo6vq/iGG2Elz3bdduBi4GjtowJ7vW5LXAi8D1pZ7H3+fBXzPafe2A6tmnK9sZZL6wOfj4PgMHQyfr0Pj81djJcm+9JIXF1fVlcN+f9fAWKCq+nxV/VBVHV5Vh9Mb5nWUyYs9S7J6xuk64G9HFcs4SXICvfVWXlFVj406Hk2sm4HVSY5Ish9wGrB5xDFJS4LPx8HxGaox4PNXYyO9v+BfCNxVVb87ihhMYGiY3pvk9iS30ZuCM/Rtd8bU7wNPAq5tW+z9wagDGhdJfi7JNuCngD9Lcs2oY+qqtsjd64Fr6C3IdHlV3THaqLovySXAXwPPTLItyVmjjkljyefj4PgMHQCfr/3j83dwfEYPxHHAa4Dj2+/UW5OcNMwA4kg6SZIkSZLUdY7AkCRJkiRJnWcCQ5IkSZIkdZ4JDEmSJEmS1HkmMCRJkiRJUueZwJAkSZIkSZ1nAkOSJEmSJHWeCQxJkiRJktR5JjAkLUiSb+3Fay9K8kCS22eVn5DkC0m2JjlnRvnjk/xFkn3a+a8k+YN2vG+SP0qyKcm+u3iv/ZLckGTZYuOVJEnDMcr+xTyub59C6ggTGJKG6SPACTMLWufhg8CJwJHA6UmObNWvA66sqn9u588GbkvyZODPga9U1fqq+u7sN6qq7wDXAf9hEF+IJEnqjI+wd/2L3bJPIXWHCQxJi5LkrUlubx9vnlH+X9tfO/4yySVJ/vPOuqq6AXho1qWOBrZW1T2tg3ApsK7VvRq4akbb57TXXw/8cVWdu4cw/7RdQ5IkjYFR9C+SvCLJx2fF8atJPjCjyD6F1AEOg5K0YEmmgDOBY4AANyX5C3q/U34eeC6wL3ALML2Hy60A7p1xvg04Jsl+wNOr6u9m1D0b+ADwuqr6X/MI9XbgJ+fRTpIkjdgI+xfvAU6f9fr/195zJ/sUUgc4AkPSYvx74E+q6h+q6lvAlcBPA8cBV1XVt6vqm8B8kgxzWQ48svMkySrgW8DngcNmN05yVJIvJ3n8zrI2NPQ7SZ60F3FIkqThGEX/4rnA46rq9iQ/muRXW9W+QO1sZ59C6gYTGJJGbTuwasb5ylb2j8ABM8qfDXyO3vDNc5M8f9Z1zgB+A3jVrPL9gW/3M2BJktR58+1fPI9/Hc3xs8DqdnwkvX7HTPYppBEzgSFpMf4PcHKSA5M8Afi5VvZXwMuTHJDkicDL5nGtm4HVSY5owzpPAzZX1cPAPkl2djKeA3y+qu4Dfgm4LMlTAFqb5cDFwMt3XjjJIcCDu1rkU5Ikdc4o+hePA57YFv18JfCkNprztcDHdl7MPoXUDSYwJC1YVd1Cb8XvvwFuAj5cVZ+tqpuBzcBt9HYJ+Tzw6M7XJbkE+GvgmUm2JTmrqnYArweuAe4CLq+qO9pLPkVvOCn0RmB8vr3/tcDlwEWt7ueBHwZ+H/ixJD/Wyl8I/Fl/v3pJkjQII+pfXA08HbgV+APgJ4AtwMYWz072KaQOSFXtuZUkzVOSJ1bVt5IcCNwAbJjVAVjItY4C3lJVr9lDuyuBV1fVPyZ5Tjt+Rys/p6q+uJj3lyRJ3TCK/sWs19inkDrAXUgk9dvGts/6AcCmxXYuoPeXmCSfSbLP7vZqr6pXzji+DbitDRf9UzsakiRNhKH3L3ayTyF1hyMwJEmSJElS57kGhiRJkiRJ6jwTGJIkSZIkqfNMYEiSJEmSpM4zgSFJkiRJkjrPBIYkSZIkSeo8ExiSJEmSJKnzTGBIkiRJkqTOM4EhSZIkSZI67/8DV8e5zv8OlF4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x667.44 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "full_ks = pd.read_csv('./data/ath.ks.tsv', sep='\\t', index_col=0)\n",
    "anchors_ks = pd.read_csv('./data/ath_anchors.ks.tsv', sep='\\t', index_col=0)\n",
    "plot_selection(dists=[full_ks, anchors_ks], title='$A. thaliana$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Mixture modeling\n",
    "\n",
    "As an illustration I will use th `wgd` library directly here, note however that the '`wgd mix`' command in the `wgd` CLI provides all tools for performing the analysis in this section without the need to code. Use for example the command below:"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "wgd mix my_ks_distribution.tsv -n 1,2 --method bgmm "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If all goes correct, you will now have the data frame with for each pair the posterior probability to come from one of the two components. Anyway, I now proceed by performing the analysis directly in the notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = filter_group_data(df=anchors_ks, min_ks=0.05, max_ks=3)\n",
    "X = get_array_for_mixture(data)\n",
    "models = fit_bgmm(X, n1=1, n2=3, n_init=10, max_iter=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4AAAADXCAYAAABLVz+FAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xdc1uX+x/HXxRJxgAMSFXPhAFdCuMhFqWmaqxyojdPR0qZlp2ycysrf6bTTMkst01xpSrnLmaAoWq40EVRw4pYtcP3+UDxoKPPmusfn+Xjcj+C+v/f3eufJz/l+vuO6lNYaIYQQQgghhBD2z8l0ACGEEEIIIYQQZUMaQCGEEEIIIYRwENIACiGEEEIIIYSDkAZQCCGEEEIIIRyENIBCCCGEEEII4SCkARRCCCGEEEIIByENoBBCCCGEEEI4CGkAhRBCCCGEEMJBSAMohBBCCCGEEA7CxXSAoqpevbquW7eu6RhCiFIUExNzWmvtbTpHSUl9EsL+2EN9ktokhP0pSW2yuQawbt26bNu2zXQMIUQpUkodNp2hNEh9EsL+2EN9ktokhP0pSW2SW0CFEEIIIYQQwkFIAyiEEEIIIYQQDsJiDaBSarpS6pRSancB292plMpSSg20VBYhhMhL6pMQwhpJbRJClAVLXgH8Buhxqw2UUs7Af4BVFswhhBA3+gapT0II6/MNUpuEEBZmsQZQa70BOFvAZk8BC4FTlsohhBA3kvokhLBGUpuEEGXB2DOASqlaQD/gi0JsO1IptU0ptS0pKcny4YQQDk3qkxDCGkltEkKUBpPLQHwM/EtrnaOUuuWGWuupwFSA4OBgXQbZrFJMTMx1vwcFBRlKIoTdk/pUBFKbhCgzUpuKKG99ktokxBUmG8BgYO7VAlYd6KmUytJaLzaYSQghQOqTEMI6SW0SQpSYsQZQa10v92el1DfAz1LAhBDWQOqTEMIaSW0SQpQGizWASqk5QGegulIqEfg34AqgtZ5iqXHthdxSJYTlSH0qPqlNQliO1KaSkfokROFYrAHUWg8pwrYPWyqHEELcSOqTEMIaSW0qHVproqKimDt3Lunp6fTs2RMfHx/TsYSwGiafARRCCCGEEKLUJCcn8/zzzxMTE4NSCq01X331FePGjZMrgkJcZWwZCCGEEEIIIUpLVlYWTz/9NDt27OD+++9n165dTJs2jVatWvH222/z4Ycfmo4ohFWQBlAIIYQQQti8adOmsXPnTgYPHswrr7xCYGAgLVu25JNPPqF169aMGzeOFStWmI4phHHSAAohhBBCCJsWFxfHt99+S6dOnXj66adxcvrfIa6LiwsTJkzA29ubhx9+mNTUVINJhTBPGkAhhBBCCGHT3nrrLZRSvPTSS7i4/H2Ki9tuu4333nuPkydP8txzzxlIKIT1kAZQCCGEEELYrLi4OGbNmkXHjh3x9va+6XbDhw+nc+fOzJgxg/j4+DJMKIR1kQZQCCGEEELYrLfffhutNeHh4bfcTinFe++9R05ODm+88UbZhBPCCskyEFagJAuX5v2uTG8shChtxa0xsiCzEMKScmvM5cuXmT9/PnfeeSfNmjUr8HtOTk707duX2bNnM2jQIHr27GnpqEJYHbkCKIQQQgghbNLKlStJSUlhxIgRhf7OiBEj0Foza9YsCyYTwnpJAyiEEEIIIWzSkiVLqFixIg899FChv1OrVi1atmzJ4sWLSU5OtmA6IayTNIBCCCGEEMLmpKSksGfPHsLCwqhQoUKRvvvAAw+QlpbGp59+aqF0QlgvaQCFEEIUS3Z2NsePHyctLc10FCGEA4qKiiIzM5NevXoV+btdu3alRo0azJgxwwLJhLBu0gAKIYQokvPnz/Paa6/RpUsXatasSefOnXn00UeJjo42HU0I4UCWLVtGpUqVaNGiRZG/6+Liwt13301sbCyRkZEWSCeE9bJYA6iUmq6UOqWU2n2Tz8OVUjuVUruUUpFKqZaWyiKEEHlJfSq+w4cPEx4ezvLlywkICGDy5MmEh4dz8OBBxowZw+TJk01HFMJmSW0qvKysLKKjo2natGm+C78XxqBBg3BxcWHSpEmlnE4I62bJK4DfAD1u8Xk80Elr3RyYAEy1YBYhhMjrG6Q+Fdnp06cZPXo058+f591332XKlCmMHj2ap59+msWLF9O8eXNmzJjBiy++aDqqELbqG6Q2FUp0dDTp6emEhoYWex++vr706tWLlStXkpGRUYrphLBuFmsAtdYbgLO3+DxSa33u6q+bgdqWyiKEEHlJfSqe8PBwTp06xeuvv063bt2u+6xKlSpMmTKFVq1a8d///pdFixYZSimE7ZLaVHjr169HKUXnzp1LtJ8nnniCs2fPsnTp0tIJJoQNsJZnAP8BLL/Zh0qpkUqpbUqpbUlJSWUYSwghpD4BREZGsmrVKgYMGED37t3z3cbNzY3PPvuMhg0b8vHHH5OYmFjGKYVwKA5dm7Zv307t2rWpWbNmifYTFhZ27QSWEI7CeAOolOrClSL2r5tto7WeqrUO1loHe3t7l104IYRDk/p0RVZWFu+//z5NmjRh7Nixt9y2fPnyrFq1CqUU77zzDlrrMkophONw9Np08eJFDh8+TPPmzUu8LxcXF0JDQ1m/fj0XL14shXRCWD+jDaBSqgXwNXC/1vqMySxCCJGX1Kf/+e677zhy5AgTJkzAzc2twO3r1avHU089xdatW5k5c2YZJBTCcUhtgg0bNpCTk0OfPn1KZX/h4eFkZmZKvRIOw1gDqJSqAywChmut/zKVQwghbiT16X+ysrL4/vvvqVWrFgMGDCj09/r160e9evX4+uuvOXfuXMFfEEIUSGrTFb/++ivlypUrlSuAAP3798fLy4u5c+eWyv6EsHaWXAZiDhAFNFZKJSql/qGUelwp9fjVTV4HqgGfK6V+V0pts1QWIYTIS+pT4UVERHDu3DmGDRuGUqrQ33N2dmbcuHGkpaXx6aefWjChEPZDalPh/Pzzz9x+++2UK1euVPbn6upKWFgYW7Zs4ezZm87BI4TdKN7CKYWgtR5SwOePAY9ZanwhhLgZqU+F98MPP+Dp6cn9999f5O+GhIQQHBzM8uXLOXDgAP7+/hZIKIT9kNpUsOTkZA4ePEiPHrdaLaPohg8fzsKFC1myZAmPPPJIqe5bCGtjfBIYIYQQ1ikmJoa//vqLbt26FerZv/y88MIL5OTkMG7cuFJOJ4RwRGvXrkVrzR133FGq++3duze+vr4sW7asVPcrhDWSBlAIIUS+FixYgFKK8PDwYu+jYcOGdOnSheXLl3P8+PFSTCeEcERr1qwBoG3btqW6XycnJ+6//36WL19OWlpaqe5bCGsjDaAQQoi/0VqzYMEC7rzzTmrXLtla00899RRZWVl8+OGHpZROCOGotmzZgq+vb4nX/8tP27ZtSUlJYeHChaW+byGsicWeART5i4mJue73oKAgQ0mEEOJ6eevT5s2biYuLY9iwYSXeb+3atXnwwQf5/PPPGTt2LL6+viXepxDCceTWJq01e/bsoVOnThYZZ8CAATzxxBMsWLCgVGqfENZKrgAKIYT4m0WLFuHq6krHjh1LZX9jxowhNTWVt956q1T2J4RwPAkJCVy8eJH77rvPIvuvWLEiISEhrFu3juzsbIuMIYQ1kAZQCCHEdbTWxMTEEBwcjJeXV6nsMzQ0lKCgIObOnUt6enqp7FMI4Vh27NgBQLt27Sw2Ru/evbl48SKrV6+22BhCmCYNoBBCiOvExMRw4cIFevXqVar7ffLJJzl//jxffvllqe5XCOEYoqOjcXV1pX79+hYbY8iQITg7OzNnzhyLjSGEadIACiGEuE7ume+hQ4eW6n6HDx9OjRo1mDJlSqnuVwjhGA4cOICvry8VKlSw2Bg1a9akQ4cObNmyxWJjCGGaNIBCCCGus3XrVvz8/KhXr16p7tfZ2ZkRI0awb98+oqOjS3XfQgj7lp2dTWJiYqnXpfwMGDCA/fv3ExcXZ/GxhDBBGkAhhBDXnD59miNHjtC9e3eL7H/8+PFUqFCBL774wiL7F0LYp4MHD5KZmUmTJk0sPla3bt0A+P777y0+lhAmSAMohBDimsjISAC6du1qkf17enoybNgw5s2bx+nTpy0yhhDC/uROANOyZUuLj9W4cWNq1KhBRESExccSwgRpAIUQQlyzatUqPD098ff3t9gYQ4YMIS0tjc8++8xiYwgh7EtCQgKurq40b97c4mMppQgNDWXHjh1cunTJ4uMJUdakARRCCAFcWf5h9+7d+Pv7o5Sy2DgdO3akbt26zJ4922JjCCHsy8GDB/H396d8+fJlMl6fPn3Iysrixx9/LJPxhChLLqYDiNIVExNz3e9BQUGGkgghbM2ePXtITk4mODjYIvvPW5+6dOnCjBkzWLNmjcVuNxVC2AetNXv37iUsLMwi+8/v2Klfv36UK1eOxYsXM2LECIuMK4QpFrsCqJSarpQ6pZTafZPPlVLqU6VUrFJqp1KqtaWyCCFEXlKf8rd+/XoAOnXqZPGxBg0ahIuLC5MnT7b4WELYCqlN+UtISCAlJYVq1aqV2ZgVK1YkKCiI3377Da11mY0rRFmw5C2g3wA9bvH5vYD/1ddIQKaEE0KUlW+Q+vQ327dvx8vLy6LP/+WqXr06nTp1YsWKFaSkpFh8PCFsxDdIbfqb33//HYDAwMAyHfeRRx4hKSmJ3bvz7ceFsFkWawC11huAs7fY5H5gpr5iM+CllPK1VB4hhMgl9envsrOzOXDgAC1atCizMceOHUtqaqo8YyPEVVKb8rd3716UUmUyAUxePXv2BODnn38u03GFsDSTk8DUAhLy/J549b2/UUqNVEptU0ptS0pKKpNw1ubPP/9k4cKFfP3118yePZsNGzaQnp5uOpYQ9srh6tOuXbtITU3l7rvvLrMx7733Xho2bMi0adPKbEwhbJzD1SaA+Ph4qlSpQtWqVct03Jo1a9KwYUOZsErYHZuYBEZrPRWYChAcHOwwN2JnZmaybNkyvv/+e+Li4v72+UsvvUTnzp0ZPHhwmZ61F0L8j73Upw0bNgDQqlWrMhtTKcUDDzzAxIkT2bt3LwEBAWU2thD2zl5qE0BiYiI1a9a06OzEN9OhQwdmzpzJqVOn8PHxKfPxhbAEk1cAjwJ+eX6vffU9AWzdupUBAwbw9ttv4+rqysSJE5k0aRILFizg22+/5eWXX6Zv375ERkby6KOPMmrUKLZv3246thD2wuHqU0REBJUrV8bXt2zvJhs6dCgAn3/+eZmOK4SNcrjalJmZyenTp42d6O7Xrx9aa+bNm2dkfCEswWQDGAGMuDqjVVvggtb6uME8VmPSpEk8+eSTJCcn8+qrrzJr1ixeeukl2rZtS7169QgMDGTAgAHMmzeP5cuXM2rUKHbt2kW7du2YMmWKzFYlRMk5XH36448/aNCgQZmfYW/WrBmBgYH8+OOPUruEKJjD1abY2Fiys7Np2rSpkfHvvfdeKlasKM8BCrtiyWUg5gBRQGOlVKJS6h9KqceVUo9f3WQZEAfEAl8Boy2VxZZMmzaNp556Cj8/P+bOnUvfvn1veUBWvnx5/vnPfzJ//nyaNWvG119/zRNPPEFycnIZphbCtkh9ut6BAwc4ffp0mU+wkGvw4MEcO3aMlStXGhlfCGshtenvNm/eDICfn18BW1qGm5sbISEhREVFkZ2dbSSDEKXNYs8Aaq2HFPC5BsZYanxrUNRF2X/++WemTJlC69at+fjjj/Hw8Cj0WLVr12bbtm0MGTKEefPmMWLECL766qti5RbC3kl9ur4+zZ07F4CQkBAjWUaNGsVbb73FV199RY8et5oBXwj7JrXp78dO0dHRANStW9dAmiv69u3LmjVriI6Opl27dsZyCFFaTN4CKvI4ePAg//d//0f9+vVZs2ZNkZq/XEopxo0bx/jx4zl58iSjRo3i+HG7vjNECFEKfv/9d9zc3Gjd2sya0t7e3nTp0oVff/2VzMxMIxmEENZp//79VKlShYoVKxrLEB4ejpOTE8uXLzeWQYjSJA2gFbh06RL/+te/8PDwYPLkyXh6epZof/379+ezzz7j5MmTtG/fnr/++quUkgoh7NHhw4dp0aIFbm5uxjI8++yzXLhwQQ6whBDXiYuLo2bNmkYzVK1alTZt2rB48WKjOYQoLdIAWoG3336bw4cP8+9//5vq1auXyj5zbyM9ceIEXbp04ehRu54kTAhRTOnp6cTFxdGtWzejOe655x5uu+02pk6dajSHEMJ6XL58mWPHjlGnTh3TUWjVqhW7du3i8OHDpqMIUWLSABq2fPlyfv31Vzp37kyHDh1Kdd9BQUHMnj2bU6dO0bVrVy5cuFCq+xdC2L7du3eTnZ1d4DPKlubi4kKXLl1YuXKl3LouhACuXP3LysqiQYMGpqMwYMAAAFkOQtgFaQANyszM5PHHH6dy5cq88sorFhmjf//+fPrppxw4cIBu3brJ8zVCiOts3boVgCZNmhhOAiNHjiQ7O5svvvjCdBQhhBXIPRlk6vnkvDp37oyXlxdLly41HUWIEpMG0KBvvvmGI0eOMHr0aLy8vCw2TkhICP/4xz+Ijo7mwQcftNg4Qgjbs3//fsqXL29sja2YmJhrr0qVKlGnTh3mz59vJIsQwrrExcUBGLkCmLc2xcTE4OzsTLt27YiOjpaT6cLmSQNoSFpaGvPmzaNJkybXbiuwpMcff5w+ffqwZMkSZs6cafHxhBC2IT4+Hj8/vzJfAD4/SikGDBjA/v37/zYVvBDC8ezcuZMaNWpQoUIF01EAuO+++0hPT5fJqoTNK1QDqJRapJTqpZSShrGUzJ49mwsXLvDll1+W2YHX+PHjCQ4OZuTIkbLgsrAbUp+KLzMzkxMnTlC/fn3TUa4ZM2YMSilZx1TYPKlNJRcfH298BtC8Bg8ejKurK7/99pvpKEKUSGGL0ufAUOCAUur/lFKNLZjJ7l24cIFZs2bRsWNHOnbsWGbjuri4MHHiRDw9PRkyZIjMZCXshdSnYsqdAMbU7Z/5adCgwbXJYHJyckzHEaIkpDaVQHZ2NidPnqRRo0amo1xTtWpVOnXqxLJly0xHEaJECtUAaq1/0VqHA62BQ8AvSqlIpdQjSilXSwa0R9OnTyc5OZkhQ4aU+dhVqlRh/vz5JCcn069fP7Kysso8gxClSepT8SUkJAAQHBxsOMn1Hn30UQ4dOiRn2YVNk9pUMvHx8WRlZREQEGA6ynU6d+7M3r172bdvn+koQhRboW9LUEpVAx4GHgN2AJ9wpaittkgyO5WWlkZERAT+/v7ceeedRjJ06tSJV199lR07dvD0008bySBEaZL6VDyxsbG4u7vTsGFD01Gu07dvX8qXL8/7779vOooQJSK1qfj+/PNP4Mr6e9YkLCwMgDlz5hhOIkTxFfYZwB+BjYAH0Ftr3UdrPU9r/RRQ0ZIB7c3s2bO5dOkSDz30kNEcr732GmFhYUyZMoVFixYZzSJESUh9Kr4//viDunXr4uzsbDrKdSpUqEDHjh355ZdfuHTpkuk4QhSL1KaSiY+PB6zvDoU2bdrg4+MjcykIm1bYK4Bfaa0DtNYTtdbHAZRS5QC01tb1N9OKaa1ZtGgRNWvWpHv37kazKKVYsGABderUYcyYMZw8edJoHiFKQOpTMWRnZxMbG4uvr6/pKPl69NFHSUtLY8aMGaajCFFcUptK4MyZM/j4+FhdjVJKERoayo4dO0hJSTEdR4hiKWwD+HY+70UV9CWlVA+l1H6lVKxS6qV8Pq+jlFqrlNqhlNqplOpZyDw2aenSpZw6dYp+/fpZxZTrVapUISIigvPnz/PQQw+htTYdSYjiKHJ9ktoEBw4cIDMzk8aNrXNeigEDBlCtWjVmzZplOooQxSXHTiUQHx9PvXr1TMfIV58+fcjMzOTHH380HUWIYrllA6iUqqGUCgLKK6XuUEq1vvrqzJVbGm71XWdgMnAvEAAMUUrd+CTvq8B8rfUdwGCuzJhlt7777jsqVarEAw88YDrKNS1atOCdd95h5cqVvPbaa6bjCFFoxa1PUpuu+OOPP4ArNcAaOTs7c9999xETE8OhQ4dMxxGi0OTYqeRycnKIjY2lWrVqpqPkq3///ri5ubF6tTzKKWxTQVcAuwPvA7WBD4EPrr7GAuML+G4IEKu1jtNaZwJzgftv2EYDla/+7AkcK3x025KUlMSiRYvo06cPFSta163/zz33HMHBwbz33ntERRV4clIIa1Hc+iS1iSsTLDg5OdG8eXPTUW7qiSeeICcnh/nz55uOIkRRyLFTCSUkJJCZmWlVawDmValSJbp168aGDRvk7ilhk27ZAGqtv9VadwEe1lp3yfPqo7UuaOaQWkBCnt8Tr76X1xvAMKVUIrAMeCq/HSmlRiqltimltiUlJRUwrHWaNm0aWVlZ9O3b13SUv1FKMW/ePNzd3Rk6dCipqammIwlRoBLUp1KrTWC79engwYP4+PhQvnx501Fuqk2bNgQHBzN37lzTUYQoNDl2Krm9e/cC4O/vbzjJzfXq1YtDhw6xZ88e01GEKLKCbgEddvXHukqpsTe+SmH8IcA3WuvaQE/gO6XU3zJpradqrYO11sHe3t6lMGzZysnJ4ddff6Vly5ZWez97/fr1+fDDDzl06BCjRo0yHUeIAlm4PhWqNoHt1qeTJ09yxx13mI5RoOHDh7Njxw62bNliOooQhSLHTiV34MABAKtbAzCvLl26ADBz5kzDSYQouoJuAa1w9Z8VgUr5vG7lKOCX5/faV9/L6x/AfACtdRTgDlQvMLWN2bBhA+fOnSM8PNx0lFt67LHH6N27N7Nnz2bNmjWm4whRkOLWJ4evTcePH+fs2bNWfXCVq2/fvjg5OTF58mTTUYQoLDl2KqG4uDgqVKhAjRo1TEe5qcaNG+Pn58eqVatMRxGiyFxu9aHW+sur/3yzGPveCvgrpepxpXgNBobesM0RIAz4RinVlCtFzLbuUyiEJUuW4OrqymOPPUZcXJzpOLc0a9Ys7rzzTh566CF27txJlSpVTEcSIl8lqE8OX5t+/fVX4MqVf2tXp04d7rzzTn7++WeysrJwcbnl/20JYZwcO5XciRMnqF27tlXMmH4rHTt2ZO7cuZw9e5aqVauajiNEoRV2Ifj3lFKVlVKuSqlflVJJeW5xyJfWOgt4ElgJ/MmVGav2KKXeUkr1ubrZ88A/lVJ/AHO4cr+8XT1Nm5qaSnR0NK1bt7aJZqpy5crMnj2bEydOWP0VSyGg6PVJahPXJnuy1iUgbjR8+HDOnTvHggULTEcRotDk2Kl4tNYcP36cli1bmo5SoPvvv5/s7GypTcLmFHYdwG5a64vAfcAhoCEwrqAvaa2Xaa0baa0baK3fufre61rriKs/79Vad9Bat9Rat9Ja29119I0bN5KRkUG/fv1MRym04OBgRo4cyfLly/noo49MxxGiIEWuT45em3bu3Em1atXw8vIyHaVQHnroITw8PGRReGFr5NipGJKSkkhJSbGJOxR69+5N+fLl+emnn0xHEaJICtsA5t5z0wtYoLW+YKE8dmfFihV4e3tfe1jYVnz00Uc0btyYV199lX379pmOI8StSH0qov3791OnTh3TMQqtYsWKhIWFsWHDBi5dumQ6jhCFJbWpGHLXKLWFWyrd3d3p0qULW7ZsITs723QcIQqtsA3gz0qpfUAQ8KtSyhtIt1ws+3D69GkiIyO55557cHZ2Nh3nOjExMde9buTm5sbcuXPJzs5m0KBBUtiENZP6VASnTp0iKSnJaqdXv1ltGjt2LBkZGfz4448G0wlRJFKbiuGvv/4CoFGjRoaT/F1+9Wn48OGcPn1a1lEWNqVQDaDW+iWgPRCstb4MpPD3hUnFDX766Seys7Np37696SjF0qpVK1555RV27tzJ+PEFrV0rhBlSn4omdzmFwMBAw0mKplOnTtSvX1+mXBc2Q2pT8cTFxeHh4WG1i8DfqGfPnri6ujJr1izTUYQotMJeAQRoAgxSSo0ABgLdLBPJfqxduxYvLy/atGljOkqxvfrqq3Tu3JlPPvmE3bt3m44jxM1IfSqk2NhYANq2bWs4SdEopejbty9r1qy5doVACBsgtamIEhISuO2223ByKsohqjmVK1cmMDCQJUuWYEdz8Qg7V9hZQL8D3gdCgTuvvoItmMvmnTp1in379tGmTRurn8b4VpRSzJs3j8qVKzNs2DDS0tJMRxLiOlKfimbHjh34+vpSrVo101GKbPjw4WitmTRpkukoQhRIalPR5c4AWqtWLdNRiqRHjx6cOHGC7du3m44iRKEUdkGlYCDAnqYZtrRvv/2WnJwcevbsaTpKifn4+DB16lT69evHqFGj5BYsYW2kPhXBhg0baNCggekYxdKqVSsCAwNZsGABn3zyiU2fXBMOQWpTEZ08eZK0tDSaNGliOkqRDBkyhP/85z/MmjWLoKAg03GEKFBhr6/vBmpYMoi9iYiIwNPT02af/7tR3759ue+++5g1a5ZMwiCsjdSnQrp06RIJCQlWOblCYY0YMYITJ05IHRK2QGpTEe3duxeA1q1bG05SNM2aNaNBgwb8/PPPpqMIUSiFbQCrA3uVUiuVUhG5L0sGs2UXL14kOjqaHj162NUZ6pkzZ+Lj48PIkSNJSkoyHUeIXFKfCmnTpk3k5OQQHGy7d6GNGjWK8uXL88UXX5iOIkRBpDYV0c6dOwFsYg3AvJycnOjatSuxsbHEx8ebjiNEgQp7C+gblgxhb5YuXUpmZibdutnXs95VqlRh2rRp9O7dm+HDh7NixQrTkYQAqU+FtnnzZgBCQ0PJzMw0nKZ4PD096devHwsXLuTcuXNUqVLFdCQhbuYN0wFszW+//Ya7uzteXl6moxTZ008/zdSpU1myZAnPPvus6ThC3FJhl4FYDxwCXK/+vBWQJ11vYsaMGXh5edG8eXPTUUpdr169ePTRR1m5ciXz5s0zHUcIqU9FsGPHDjw8PGjWrJnpKCUybtw4MjIyZNp1YdWkNhXdgQMH8PHxwcWlsNcnrEdgYCDNmjWT29OFTSjsLKD/BH4Avrz6Vi1gsaVC2bKsrCyioqIICQmxmSmMi2ry5MkphiKkAAAgAElEQVQEBwczevRojh8/bjqOcHBSnwovLi6Oxo0b2/yt6a1ataJVq1Z89tlnMu26sFpSm4ru0KFDNjcDaF6dOnVi48aNcmwkrF5hO5QxQAfgIoDW+gDgY6lQtmzVqlUkJyfbxeyfN1OuXDlmzZpFWloaDzzwgByACdOkPhVCVlYWsbGxdO7c2XSUUtGnTx8OHDjAL7/8YjqKEDcjtakITp48ycWLF6lbt67pKMXWq1cvtNZ88803pqMIcUuFbQAztNbXHhhRSrkActSfj0WLFqGU4sEHHzQdxaIaN27MCy+8wKZNm3j99ddNxxGOTepTIezZs4f09HSbm13vZp588knc3Nz47LPPTEcR4makNhXBtm3bAGjYsKHhJMUXFhZGrVq1mDt3rukoQtxSYRvA9Uqp8UB5pdQ9wALgp4K+pJTqoZTar5SKVUq9dJNtHlRK7VVK7VFKfV/46NZp3bp1NGnSBF9fX9NRSkVMTMx1r7zefPNNgoKCeO+999i6dauhhEIUvT45Ym1atWoVAP7+/oaTlI4jR44QFBTEypUrWbt2rek4QuRHjp2KIC4uDoAWLVoYTlJ8bm5udO3alZ07d7J06dK/HTcJYS0K2wC+BCQBu4BRwDLg1Vt9QSnlDEwG7gUCgCFKqYAbtvEHXgY6aK0DAZueNunkyZMcPHiQ7t27m45SJpRSzJs3D1dXV4YMGUJGRobpSMIxFak+OWJtgisnc1xdXbnjjjtMRyk1AwcOJDMzUyakEtZKjp2K4MCBA1SqVMmmbwEFGDp0KICsCSisWmFnAc3hyoPLo7XWA7XWX+mCH/wKAWK11nFXb4GYC9x/wzb/BCZrrc9dHedU0eJbl9xlEUaMGGE4Sdlp0KAB7733HgcPHuTJJ580HUc4oGLUJ4erTXBlgeW6devi5uZmOkqp6dixIzVr1mTZsmXyLLKwOnLsVDR//PGHXUxS1blzZ/z8/OTOBGHVbtkAqiveUEqdBvYD+5VSSUqpwjz0VQtIyPN74tX38moENFJKbVJKbVZK9bhJjpFKqW1KqW3WvAD54sWL8fX1pVWrVqajlKnRo0fTu3dvpk+fTlRUlOk4wkGUoD6VWm26msPq65PWmoMHD9K0aVPTUUqVUopHH32UhIQENm7caDqOEIAcOxXXrl277GJdT3d3d/r27cv+/ftJTEw0HUeIfBV0BfA5rsxgdafWuqrWuirQBuiglHquFMZ3AfyBzsAQ4Cul1N9W/9RaT9VaB2utg729vUth2NKXmZnJqlWraNWqlc2fvSqO7777Dj8/P4YPH87FixdNxxGOwZL1qVC1CWyjPu3evZvU1FS7uv0zV48ePahcuTIff/yx6ShC5JJjpyI6e/Ys586do3HjxqajlIp77rkHgGXLlhlOIkT+CmoAhwNDtNbxuW9oreOAYUBB9zkeBfzy/F776nt5JQIRWuvLV8f4iytFzeYsX76c1NRUevXqZTqKEZ6ennz77bfExcUxZMgQ03GEYyhufXKo2gRXGkC4csukvXF3d+euu+4iIiKCgwcPmo4jBMixU5HlTpbSrFkzw0lKR/Xq1alVqxbLly83HUWIfBXUALpqrU/f+KbWOglwLeC7WwF/pVQ9pZQbMBiIuGGbxVw5g4VSqjpXbmuIK0Ruq7N48WKcnJzsfvmHW+nUqRMPPfQQy5Yt4/333zcdR9i/4tYnh6pNcGUJCGdnZ9q1a2c6ikUMHjyYnJwc3nvvPdNRhAA5diqy7du3A9jNMjVubm60a9eOhIQEduzYYTqOEH9TUAOYWczP0FpnAU8CK4E/gfla6z1KqbeUUn2ubrYSOKOU2gusBcZprc8ULrp1WbduHU2bNsXab7MoDTdbFgLgyy+/JCAggFdffZXo6GgD6YQDKVZ9crTaBLBlyxaaNGlC+fLlTUexiKZNmxIUFMT3339PVFTUTeuTEGVEjp2KaNeuXbi5udnNFUCA3r174+TkxBdffHHLJbWEMKGgBrClUupiPq9LQPOCdq61Xqa1bqS1bqC1fufqe69rrSOu/qy11mO11gFa6+Zaa5tcOTM2NpZDhw4RFhZmOopxbm5uLF68GDc3NwYOHCjPAwpLKnZ9cpTalGvbtm3UrFnTdAyLGjNmDMnJySxcuNB0FCHk2KmIjh8/bncnqRo3bkyjRo2YP38+ly9fNh1HiOvcsgHUWjtrrSvn86qktS7oNgaHsW7dOuB/a784On9/fz7//HMSExMZM2aM6TjCTkl9KpxDhw5x/vx5u5+deNiwYVSrVo358+fLkhDCKKlNRbdv3z67m6TKxcWFrl27cuHCBTZt2mQ6jhDXKexC8OIWVqxYQa1atQgJCTEdxWoMGzaMl156iVmzZvHdd9+ZjiOEw8pdHsHe65OLiwvDhw8nISFBbrESwoacO3eOY8eOUb9+fdNRSt3AgQO57bbbZFF4YXVcTAewdRkZGSxfvpyBAwfa/PIPNx40BQUFlWh/b731FpGRkYwcORJ/f3/atm1bov0JIYou91lcW54BtLC16YEHHmDmzJnMnDmT4ODgsogmhCihrVu3AuDj42M4SfHcqj5VrlyZYcOG8fHHH3Pu3Dm7WOdQ2Ae5AlhCP//8M6mpqXTo0EEe8r2Bi4sL06ZNw9XVlQEDBnD69N8mRRNCWNjOnTupUqXKtStj9lybypUrR+/evYmMjLy29IUQwnrFxMQQEXFlklMPDw/DaSzj3nvvJTs7W64CCqsiDWAJLV68GGdnZwYMGGA6ilVq0KAB06dP58SJE/Tu3Zvs7GzTkYRwKImJiTRt2tR0jDIzcOBAXF1d+eqrr0xHEUIUQmxsLG5ubjRo0MB0FIto3rw5vr6+LFy4UJ5PFlZDGsASWr9+PYGBgVSrVs10FKs1cOBAXnzxRTZv3szo0aNNxxHCYVy4cIG4uDhatmxpOkqZ8fX1pUOHDkRFRXHo0CHTcYQQBTh8+DA+Pj64u7ubjmIRPj4+dOjQgcTExGvrHQphmjSAJbBv3z4SEhK4++67TUexeu+++y7dunVj6tSpzJ8/33QcIRxCZGQkAA0bNjScpGw9+uijaK157bXXTEcRQhQgMTGRWrVqmY5hUb1798bd3Z05c+aYjiIEIA1gieTetz5w4EDDSayfUooffviBli1b8thjj7Fv3z7TkYSwe+vXrwew21urbiYgIICgoCDmz5/P0aNHTccRQtzE+fPnSU5OtvtJmxo2bEjLli3ZuHGjzIcgrII0gCUQGRnJ7bffLrNbFlKlSpX46aefcHd3p1evXpw4ccJ0JCHs2o4dO/D09KR27dqmo5S50aNHk5WVxUcffWQ6ihDiJg4ePAhAkyZNDCexrHLlyvHAAw+QnZ3N4sWLTccRQhrA4kpPT+eXX36hZ8+eNr/8Q1ny8/Nj5syZHD58mO7du5Oenm46khB2688//6Rx48amYxjRokULhg4dyueff86pU6dMxxFC5CN3tl57vwUUoHPnzrRt25ZFixaRlZVlOo5wcNIAFtPixYtJSUkhNDTUdBSb06NHD95991127txJ//79ZVYsISzg4sWLJCYm0rx5c9NRjHnllVfIyMhg/PjxpqMIIfKxf/9+ypUrh7e3t+koZaJfv36cOnWKBQsWmI4iHJw0gMUUERGBi4sLvXr1Mh3FJr344os89thjLF++nKefftp0HCHszm+//YbWmpCQENNRjGnSpAl33XUXs2bN4vjx46bjCCFuEB8fj4+PD+XLlzcdpUzUrFmTatWqMXHiRDn5LYySBrCY1q9fT4sWLfD09DQdxbjcxaWLusj01KlTCQsLY9KkSUybNs2CCYVwPPHx8QDcc889hpOYFR4eTmZmJo899liR6pMQwrK01hw9etQhbv/MVb9+fdq1a8euXbv4/PPPi3XsJERpsGgDqJTqoZTar5SKVUq9dIvtBiiltFLKJqaB2rVrF8eOHSMsLMx0FJumlCIiIoI777yTMWPGsGHDBtORhIOw19qU1/bt26levTp169Y1HcWo1q1b06ZNG1auXMnhw4dNxxGiQI5QnwCOHz9Oamoq9erVMx2lzJQrV47u3bvj6enJt99+azqOcGAWawCVUs7AZOBeIAAYopQKyGe7SsAzwBZLZSltuevYDRo0yHAS2+fh4cHy5cupV68e9913H2vWrDEdSdg5e65NeW3cuJGAgACZpAp46qmnUErx/vvvm44ixC05Sn0CiI6OBux/BtAbNW3alK5du7J161ZZEksYY8krgCFArNY6TmudCcwF7s9nuwnAfwCbmQ5y8+bN1K5dm9atW5uOYheqVavGypUrcXNzo1+/fmzfvt10JGHf7LY25UpJSSEuLs6hzqzfSuPGjenZsydRUVFs3brVdBwhbsXu61Ou3ObH0SbT8/LyYuDAgbi7u8vjL8IYSzaAtYCEPL8nXn3vGqVUa8BPa730VjtSSo1USm1TSm1LSkoq/aRFkJ6ezqZNm+jbt6+cWS9FderUYcWKFWit6d69O3/99ZfpSMJ+lVpturqt1dSnXJs2bSI7O9uhJ4C50dixY6lSpQrjxo2TyReENbPLY6f8/PHHH/j6+jrkXAqDBg1i+PDhrF27Vq4CCiOMTQKjlHICPgSeL2hbrfVUrXWw1jrY9FTBy5YtIy0tjZ49exrNYe2KMzFMcHAwP/74I8nJyYSFhXH06FELpxTi74pSm8C66lOuyMhIAO666y7DSaxHxYoVeeSRR1i/fj0TJkyQiReETbLVY6f8REVF4evrazqGEc7OzgwdOpSKFSvy5ZdfXnu/uJPqCVFUlmwAjwJ+eX6vffW9XJWAZsA6pdQhoC0QYe0PM8+dOxdXV1eHu2WhrISFhTF79mxOnDhB165dOXPmjOlIwv7YZW3Ka/v27VSoUIFmzZqZjmJVevfujZ+fHx999BGpqamm4wiRH7uvTwCpqakcOXKEGjVqmI5iTGpqKsHBwWzcuJHdu3ebjiMcjCUbwK2Av1KqnlLKDRgMROR+qLW+oLWurrWuq7WuC2wG+mitt1kwU4lt3LiRli1bUqlSJdNR7Fb//v2ZPXs2hw8f5u6775YmUJQ2u6xNee3fvx9/f3+5Tf0GlSpV4t///jfnz5/n008/NR1HiPzYfX0C2LFjB1pr/P39TUcxplq1arRr144KFSowefJkuTVdlCmLNYBa6yzgSWAl8CcwX2u9Ryn1llKqj6XGtaSYmBhOnDhBt27dTEexew8++CBLlizhzz//pE2bNnI7qCg19lib8kpPTyc+Pl7q1E089NBDtG7dmkWLFnHgwAHTcYS4jr3Xp1y5t6k3b97ccBJzXFxcCAwMpH379mzdupWIiIiCvyREKbHoM4Ba62Va60Za6wZa63euvve61vpv/5VrrTtb+xmsefPmATB48GDDSRxD9+7dmT59OocPHyY0NFSaQFFq7K025bVjxw4uX75M27ZtTUexSi4uLjzzzDO4ubnx5ptvkpOTYzqSENex5/qUKyYmBnd3d4e+AgjQoEEDQkJCqFGjBs8//zyZmZmmIwkHYWwSGFu0atUqateu7dBnrMra0KFDmTFjBomJibRp00ZmBxWiALlnkYOCggwnsV6BgYGMGDGCffv2MXnyZNNxhHA4f/31F/Xr16dChQqmoxjl6upKo0aN6Nq1KwcPHmTOnDmmIwkHIQ1gIV26dIm9e/fSr18/01EczrBhw5g3bx5nzpyhffv2bNtmcyc7hSgzmzdvplq1atSpU8d0FKv2z3/+kw4dOvCvf/2L2NhY03GEcBg5OTkcOHCAsLAw01GsQqNGjXj88cfp06cP06ZN48SJE6YjCQcgDWAhrV69msuXLzNgwADTURxS//79Wbp0KdnZ2dx7770yPbIQN7Fnzx4CAwNNx7B6SileeOEFnJycGDp0qNwKKkQZ2b9/P8nJydxxxx2mo1gFV1dXPDw8+OSTT9Ba8+6778qEMMLipAEspJkzZ1KpUiXat29vOorD6tq1K9HR0VSoUIEuXbowd+5c05GEsCqHDh0iKSmJ4GCbmhHemFq1avHII4+wdetW3n77bdNxhHAIv/zyCwD169c3nMS6xMfHc++99xIZGcnSpUtNxxF2ThrAQsjOzmbdunUEBQXh6upqOo5D8/f3Z9OmTdx2222Eh4fz+uuvm44khNXIPbDq1KmT4SS2wcnJiddee42QkBDefPNNNmzYYDqSEHYvMjISd3d3WrZsaTqKValevTqNGzcmICCADz/8kKSkJNORhB2TBrAQ1qxZw4ULF+jZs6fpKIIrZ+23bNlCcHAwEyZM4MEHHyQrK8t0LCGM++OPP3B2dqZLly6mo9gMHx8fJk6cSPXq1XnwwQc5deqU6UhC2LXff/+devXq4eXlZTqKVQkICKBixYr06NGDzMxMXnnlFbKzs03HEnZKGsBC+OGHH1BKyfIPZSQmJua6V36qVq3Kpk2bGDRoEAsWLKB9+/Zytkw4vL1799KqVSsqVapkOopN6dSpE+PGjePs2bMMGTKEy5cv57tdYWqTEOLmMjIyiI2NleeU8+Hs7Ezz5s1xd3fnkUceYfv27bz11luF/r7UJlEU0gAWwq+//kqjRo3w8/MzHUXk4eLiwty5c3n99deJiYmhbdu27Ny503QsIYzIyMhgy5Ytsv5fMTg7OzNq1Ci++OIL1qxZw5NPPimTMAhhAdu3bycrK0ueU76JmjVr4uPjQ506dejVqxcTJkxg5cqVpmMJO+RiOoC1S0hI4ODBgzz55JOmo9i8G89KldY6ZW+++SadO3cmPDyctm3b8t5778n/XsLhrF27lpSUFDmzXgy5talVq1aMGDGCqVOn0rRpU5599lnDyYSwL1u2bAHgnnvuMZzEegUFBeHi4kJoaCgJCQkMGDCAadOm0bBhw2ufC1FScgWwALmLKo8ePdpwEnErXbp0Yfv27TRv3pynnnqKfv36kZKSYjqWEGVm9erVAHTv3t1wEtvWqVMn7rjjDsaOHcv8+fNNxxHCrmzfvh1fX19ZAuIWypcvj6urK+XKlWPOnDl4eHjwzDPPcPr0adPRhB2RBrAAP/zwA02bNqVp06amo4gC1KhRg40bNzJ06FAWL15MQEAAv/32m+lYQpSJ3377jZo1a8rU6iUUEBDAs88+i7+/P+Hh4fz000+mIwlhN3777TdCQkJQSpmOYvWio6PZuXMn//3vf7l48SJjxozh7NmzpmMJOyEN4C0kJiayYcMG7rrrLtNRRCG5ubkxe/ZsZsyYwfnz5+nSpQvjx4+XRZ6FXcvKymLXrl20bt3adBSb5+rqSp8+fXjmmWeoW7cuAwcOlGdwhCgFR48eJT4+npo1a5qOYhMaNGjAxYsX0Vrz4YcfkpiYyOOPP87JkydNRxN2QBrAW5g5cyY5OTk8+OCDpqOIInr44YfZuXMnQUFBTJw4kR49enDkyBHTsYSwiE2bNpGWlibr/5USLy8vevXqxejRo6lTpw69e/dm7ty5pmMJYdNWrVoFIMvUFJK3tzdNmzYlNjaWevXq8cknn3Ds2DE6duzIgQMHTMcTNk4awFtYsmQJ1atXp2vXrqajiGK4/fbbiYqK4osvvmDTpk00bdqU8ePHy5qBwu7kTmLSu3dvw0nsx+23306PHj345ZdfaNeuHUOGDGHOnDmmYwlhs3755Rfc3d3p2LGj6Sg2o02bNnh4eLBt2zZatGjBpEmTOHPmDG3atGHdunWm4wkbZtEGUCnVQym1XykVq5R6KZ/Pxyql9iqldiqlflVK3W7JPEVx9uxZtm/fTlhYmNyrbsOUUjz++OPs3r2bJk2aMHHiRJo1a0ZkZKTpaMIgW65N+dmwYQP16tWjcePGpqPYlaZNm3L77bezcuVKevbsyQcffMCECRNIT083HU3YMXurT7k2b95Mw4YN8fHxMR3FZri5uRESEoKbmxuZmZm0atWKLVu2UKNGDe6++27efPNNOaktisViy0AopZyBycA9QCKwVSkVobXem2ezHUCw1jpVKfUE8B4wyFKZimL69OlkZWUxdOjQfD+31JIGjqas/hzr1avHtm3beP/993nzzTfp2LEj4eHhfPTRR1StWtUiYwrrZOu16UYZGRn88ssvDBkyBJDaVFry/jnGx8fTu3dvPDw8+OGHH9i3bx//+c9/5M9WlDp7q0+5zpw5Q1xcHL169WL79u2m49iUKlWq0KVLl2sXIxo0aMAXX3zBf/7zH9544w1++OEHZs2aZTilsDWWvAIYAsRqreO01pnAXOD+vBtorddqrVOv/roZqG3BPEWycuVKatWqJbdU2RGlFOPGjePPP//k7rvvZubMmQQEBPDtt9/KJDGOxaZr042WLVtGSkqKLKxsQXXq1MHPz49GjRrx8ssvc+zYMQYPHszHH39Mdna26XjCvthVfcq1du1aAMLCwgwnsU1KKbKysti2bRvHjh2jYsWKTJgwgXfffZeEhATuuOMO3nnnHZKSkkxHFTbCkg1gLSAhz++JV9+7mX8Ay/P7QCk1Uim1TSm1rSz+4z558iRr1qzhoYcekts/7ZCfnx8rVqxg3bp11K1bl4cffpjmzZuzePFi09FE2Si12gRlX59uFBERgZOTE/379y/zsR2Fs7Mz3bt3x8fHB6UUH3zwAcHBwTz33HN06NCBTZs2mY4o7IfNHjvdyvLly/Hy8qJ9+/ZGc9gyrTXnzp1j9erVJCcnA9CtWzd+/PFHnnnmGSIiIujduzdvvvkmO3bsQGttOLGwZlYxCYxSahgQDPw3v8+11lO11sFa62Bvb2+L5/n888/Jycm56e2fwnbExMRc98qrU6dOREZGMn36dI4fP06/fv0ICgq6tqB2WWYT1qmg2gRlX59utHHjRpo0aYKJsR2Ji4sLbdu2xcfHhyNHjjBx4kRmzZrF4cOHCQ0NpW/fvkW6te1Wtck0a84m/sfajp1uRmvN0qVLCQ0NxcXFYk8e2T1XV1fatWsHXHmeMjMzE4DKlSvz0UcfsWjRIvr378/q1atp3bo1gYGBTJgwgejo6CI/J2itNcBac9kiSzaARwG/PL/XvvredZRSdwOvAH201hkWzFNo33//PfXq1SMwMNB0FHETpVUEnJyceOSRRzh06BBjx45l//79dOvWjdDQUKKjo0sxsbAiNlubbnTo0CHi4uJkVr0yktsEhoaGUr58ecLDw4mNjeWdd95h7dq1BAUF0bFjR/7v//6PTZs2yUGKKA67qU+5du3axcmTJ2WSqlJQsWJF7rnnHpKTk4mMjLyusatVqxYvvvgiS5cuZcqUKVSvXp3XX3+dNm3aUK1aNe655x7Cw8P597//zZQpU1i4cCGXLl0q8ZVCrTVZWVlkZGSQkpLCxYsXOX/+POnp6XIV0opZ8lTMVsBfKVWPK8VrMHDdJTWl1B3Al0APrfUpC2YptOjoaGJjYxk7dqzpKKIMVa5cmQ8++ICXX36ZV155hZkzZ9KmTRv69OnD2LFjZX01+2KTtSk/s2fPRmstdyuUIRcXF6pXrw5AbGwsR44c4cUXX2T06NFMnz6dzz77jJdffhkPDw86duxI+/btqVGjBrVq3eouPiGusZv6lOvHH38E4L777jOcxD7UrFmTkJAQtm/fTnJyMl5eXtd97unpyahRoxg1ahQnT55k3bp1rFmzhpiYGDZu3EhGxvXnC5ycnPDw8KBChQqUL18e4FrzlvtSSpGdnU12djZZWVnX/fNWcygopShfvjweHh54eHhQpUoVvL29r3vVqVOH22+/nbp161KrVi25SlxGLPanrLXOUko9CawEnIHpWus9Sqm3gG1a6wiu3LZQEVhw9Vm7I1rrPpbKVBgffvghzs7OPPXUUyZjCEOqV6/Ol19+ydtvv83kyZP57LPPiIiIICAggBdeeIERI0bg7OxsOqYoAVutTfnZuHEj9evXJzQ01HQUh5SamkpsbCxnz56lW7dujB07lmeeeYapU6eyevVq1qxZw4oVK3j99depX78+LVq0IDAwkKZNm3Lp0iWqVq1K1apVKV++PNnZ2ddqi9aay5cvk5GRQUZGBmlpaaSmppKamkpKSkqBP+f3WVZWFu7u7te9PD098fb2pnr16tcdjGVkZFCuXDnDf7qOyZ7qU64lS5ZQu3ZtgoOD2b9/v+k4dqFmzZp4e3vj6uoKcNMm7LbbbmPQoEEMGnRlktjo6GiOHj3KqVOnSEpKws3NjYsXL15XO5RSnDt3DqXUtVeNGjVwdnbGxcUFZ2fnfH92cnLC2dkZV1dXsrKyOHnyJJcuXSIlJYW0tDTS0tLIzs7m0qVLHDx4kFOnTpGSknJdXmdn52sTbgUGBhIYGEhAQAABAQF4enpa9g/VwVi0zdZaLwOW3fDe63l+vtuS4xfVpUuX+PnnnwkNDaVu3bqm44gylPc2raCgILy9vXnjjTd4/vnneffdd/nqq6949NFHGT9+POHh4Tz//PP4+vrKlPs2ytZqU37OnTvHr7/+ytixY2WyKkNatGhBlSpVWLNmDYsWLaJr167UqVOHkJAQQkJCeOmllzhw4ACnTp1i06ZN7Nmzh59++ummM4fmHkzlPttTVE5OTlSoUOHa2XYnJ6drzV6VKlXIyMggKSmJ9PR00tLSuHDhAmfOnMn3Ni1vb298fX2pXbs2HTt2pFmzZgQGBlKnTh2cnKxi+gC7ZQ/1KdfRo0fZsWMHgwcPpmLFiqbj2JXc5u+vv/7i6NGj+Pj4XHvvZpydnalTpw516tQBbn7MUtCxTVZWFsuWLePChQtcuHCB5ORkvL29adeuHQEBAZw9e5Yffvjh2pjOzs44OTkRGhpK/fr1WblyJVFRUWRkZHD+/HnKlSvHmTNn8PT05OzZs+zatYu1a9dy+fLla2PWqFEDf3//65pCUXxynTWP2bNnk5KSwpgxY0xHEVaiUqVKTJw4kTfeeIPPP/+cr7/+mg8++IBJkyYxZMgQOnXqRLNmzeQAXJS53ItSUWAAACAASURBVLVK5bYqs/z8/OjXrx+rV69mxYoVDBw48Npnzs7ONGnShPDwcJ577jngyrqN8fHxrF+/njNnznDu3DnS09Px9vYmIyOD7OxsypUrh5ubG+XKlePkyZPXXbVr0aLFtdu1chu93J/d3Nyuq0VRUVGcO3eOtLQ06tatS0pKCpmZmQQHB+Pl5cXhw4fZvHkzycnJJCcnc+HCBc6dO8f58+c5c+YMhw8fJioqimXL/teLVKxYkWbNmhEcHHyt0fX395emUORr/vz5ANf9vRCly8PDg1OnTrF//37atWt37TbO0pKTk8PJkyc5fvw4lStXpn79+uTk5LBlyxaUUlSqVIlKlSoRGBhItWrVAPDy8mLo0KG4u7vne0tn9erVry2zprWmefPmpKWlXatjZ8+eZe/evcTGxvLnn39y8OBBEhISOHToEBs3bry2H19fX1q2bEmHDh0ICwsjKCgINze3Uv33t1fSAOYxc+ZMGjVqJIXKCt14ha6gbW61XVHHy93Xc889x3PPPcfGjRuZNWsWs2fP5ptvvsHb25suXbrwwAMPEBQUlO93i5utJP9OcnXSvs2ZM4dq1arJtOqG5f498/Pzw8XFhapVqxIfH096ejru7u5/2w6u/F288danotS1mJiYa7dSXbp0iZYtW3L69GlOnTrFmTNnCAwMxMfHh7NnzxIZGQlcWYj7xIkTuLm5obWmatWq1KxZk9OnTwNcy966dWsqV65MhQoViIuL4/fff8fb25tjx45de2VmZjJjxgwmTZoEQIUKFQgICLg2gVb79u2pXLlygX9mBf27C9u3cOFCmjdvzv3331/wxqJYateuTevWrfn6669Zu3Ytbdu2BQr/9+xm2x06dIijR49y+vRpate+stRkkyZNqF+/Pm5ubnTt2pWKFStea/Bya1NiYmKBY+allKJcuXLs3r37uvdzJ9tq3rz5tbH37t1LamoqGzduJCoqiqSkJLZt28aKFSt47bXXcHd3584778TPz49WrVoRFBSEp6dnsY+d7Lk2SQN41erVq4mKiuKTTz6RqzkiX7lFwcPDg5EjR/Lf//6X8ePHs3z5cubPn8/8+fNp3rz5/7d359FRlffjx9/PTDJZyE4C2XcUAgGSKKIQoigCaqV1aesfbq302+38fr9TT4+12Nbft9bq13N6tPXnwkFPrVbr0h5NlWqViJhACEsSIAmB7Ps+2ckymef3xyS3BFEWk8wM83mdM+dMMjc3n7lz72ee57nPwtq1a7n99tu/MDBbiJly8uRJDh8+zH333SdjUl2E2Ww2Ckm9vb3s3r2b+Pj4WZ1Nenh4mKKiIoqLi6mtrQXAz8+PgYEBbr31VsLCwli3bh2+vr6sXbuWkpKSaX8fFRVlFBannN5glZycTGJiIosXL2ZoaIjBwUH6+/tZuXIlWmt+97vfUVRUREtLi7E0ht1ux2QykZGRQXZ2Njk5Oaxdu9aYOEd4jhMnTlBQUMDjjz8uE3vMsri4OHJycti3bx+ff/75Rc0M3dPTQ29vr1HpaWxsZHh4mISEBNavX09kZOS0u4vnW8aZqQafefPmYTKZCAgIYPPmzWzevJmMjAysVivl5eUUFBTQ2dlJfn4+f/vb33j99dcBR55bu3Yt69evZ+PGjSQmJko5H6kAGn71q18REBDA9773PWeHItxEUFAQ999/P/fffz9VVVW8/fbbFBQU8Pzzz/Piiy+SmZnJypUrefDBB2f0/0rrufjjH/+I1lomq3JRAQEBpKSkUFVVRUtLC6GhoUbF6GL19/fT3t5OW1sbPT09pKWl4ePjg9lsZvny5URGRhIaGjrtrqPFYjEqXhfbUGAymYwuXme66aabuO6664xp39vb2zlx4gQTExMcPnyY5557jqeffhqA1NRUcnJyuOGGGwgJCZF1Kz3As88+i9lsZvPmzc4OxSMEBwdz3XXX0dHRYUziZLPZvrLy3dPTw7Fjx2hqamJ4eBiTycQ3v/lNAK6++mrjb5OSkmb/DVwEk8nE/Pnzyc7OJjs72/j9+++/bzSOVVZW8v777/P2228DjvHNl112GZs3b2bjxo1kZGR4ZEOqVACBvXv3sn//fn7wgx/IIGVxUVJTU3n44YfJzMzklVdeIT8/n127drFjxw527NhBbGyssUbY1EKuZ7oUK3aX4ntyNpvNxltvvUVaWhoZGRnODkechZeXF+np6cTHx1NaWsrnn3+O1Wpl/fr1F1wJPH78OA0NDQwODgKwcuVKY4yL2WwmOzv7rN3MZ9rZrmWllDEOMTIyksWLF09bMkdrzaeffsr+/fuprKzkr3/9Ky+99BLguGORnp7OqlWrjNlH55rkp9kzOjrKX/7yF1asWEFqaqqzw/EYPj4+xMU5lpFsa2ujuLiY5cuXn3UZmrq6Og4dOkRDQwMREREsWbKE6OhoYyIZd7lre7brOCoqii1bthhdj1NSUnjzzTepra2lpKSEffv2UVBQwCOPPIK/vz/Lly8nJyeHzZs34+XlNa0hzRnmIje5x6c7y37+85/j6+vLr3/963NvLNzebF5YSinS09NJT0/nRz/6EWVlZRw8eJDc3Fxyc3N577338Pf3Z9OmTSxatIj09HSnLo4rBSD38+6779LR0cGjjz7q7FDEOQQHB7Nu3TrCw8MpLCw0Kn91dXVnvSN48OBBuru76e7u5vLLLycrK4uBgQH8/f1JTU0lMjKS7Oxst1lc/oorruCKK65Aa01PTw/Nzc1UV1ezZ88eDh06RF5eHjt37uTRRx8lKiqKNWvWcMMNN3DDDTeQnJzM4cOHp+1P8pP7ePnll+nr6+O2226ThnUn8fX1xdvbm/379xMWFsb4+DgFBQWkpKQQFxfHwoULWbJkCVar9bwqPO48Ni4kJMTIR3feeSejo6MopYzhX0ePHuXJJ5/kySefxMvLi/j4eJYtW8add97Jxo0bv9Bj4VIoO3l8BTA3N5e9e/fy4x//+EsX6r0UPmjhHEuXLuWee+7h3nvvpauri7y8PMrLyzlw4AD/+Mc/AEd3sczMTJKSkli9ejUZGRmz1vLmzglcODzzzDMkJiaydetWyU1uIiEhwZhspbu7m0OHDtHR0UFKSgoJCQkMDAxQUFBAXl4eIyMjmEwmEhISALjyyiudGfqMUEoxf/585s+fz/Lly4mPj+euu+6iubmZI0eO0NjYSEFBAR988IExdXxMTAxJSUmsWrWKa665xjgewvWNjY3x2GOPERMTw4oVK4w8JflpbgUHB5OSkkJRURG7d++mtLSU0NBQFi1aBDjGC6elpblNg9JM8vHxISsry+iRNTExYUwo89prr7Fv3z4++OADcnNzAcd6illZWVx11VXcfffdaK3dfhyhR1cAJyYm2LZtGyEhITz++OPODkdc4sLDw/n2t78NQGZmJq+99hr5+fmUlJRQVlbGnj17eOWVV4wZ9a688kr8/f3JzMxk0aJFbp9sxNf35ptvkp+fzzPPPOM23XPEdGFhYVx++eX8/e9/59ixY+Tn52O1WsnOziYsLIyYmBgiIyPPuZ6Xuzt9PbKsrCz279+P1Wqlo6ODiooKysrKyMvLIz8/nz/84Q8EBgayevVqNm/ezIYNG0hLS5OlJ1zU9u3baWlp4b777iMyMtLZ4XgUm81Gf38/4Gh4qa6uJiQkhC1bthAVFcWpU6eIiYnBZrNNW2PP05nNZiIjI7nllluIiooCHOMju7q6eOeddygvL+ezzz5j586d/OY3vyEwMJDk5GQuv/xyli9fTnR0tPF37sKjSxDPPvssx44dY8eOHQQHBzs7HOFBlFLTFjLNysrinXfeYd++fTQ0NFBeXs727dux2WyAoytHfHw8q1evJjk5mZSUFBITE40ZuWZ7CYyZJHetLs7ExAQPP/ww4eHhPPDAA84OR1yg07tBWiwWYmJiWLlyJVVVVZw8eZL777+fsrIyWltb6e3tNdbTcnfne717eXkRERFBREQES5cu5fHHH2f37t2UlZVRVFTE0aNHOXjwIB9//DHwn6UnsrKyWLt2LTk5OcYsrHMRrzi7lpYWtm3bRnZ2Nrfffrs0XM6BiYkJ6urqKCoqorW1FaUU1113HeBYSmGqe+fpY4WbmpooKSlhbGwMm812yeSbC/VV13tYWBgbNmzgsssuAxzrIdbU1FBdXc27775LbW0tpaWlvPXWWzzyyCMsXLiQtLQ0kpKSSElJ4Z577pmRnHQ+sV4Mj60AFhcX88tf/pKbb75ZZv4ULiEpKWnaTFtLlizhlVdeoaSkhIqKCmpqanj99deNSqFSyujOER4eblQMFy9ejNZ6RmObybUNxcX57W9/S21tLU899RT+/v7ODkdcgH379pGXl2fMspednY3NZjPWp8rKyjIKaeXl5fT19WGxWBgaGqK3t5eIiAinT0ow15RSBAUFcfXVVxvdtNLS0vj3v/9NcXEx+fn5nDx5ku3bt/PCCy8AjpmZr732WmJjY0lKSmLx4sVERkZyxRVXOPOteAy73c7WrVsZHR01xgCK2dXU1ERxcTHR0dF0dHSQkJBAbGysMavll+WNBQsWkJKSQkNDAydPnmTevHkkJCQ4dU4CV2cymUhNTeU73/kOGzduBKCrq4vDhw/T2dlJWVkZFRUVfPrppwBs27aNefPmGd3Zr7zySm655RZGR0eZN2+esV9nlZ08sgJotVrZsmULJpOJZ599VlqohEvy9/dn1apVrFq1yvjdihUreO+996iurubgwYNUV1fT0tLCgQMHsNvtxnbz5s0jKCiIiIgIoqOjiY2NZc2aNWitSUxMPOuU7sJ1HT16lCeeeGJWlhURM2tkZITGxkaGhoZYvHgx4FgXMDg4mLS0NGPtvS8bd5OTk0N7ezstLS3U1dVx8uRJ4uLijLGANTU1DA0NTStAeAo/Pz9iY2OJjY3lG9/4BhkZGXz00UeUlpZy5MgRqqqqOH78uDFuBxwF4EWLFrF48WIyMjLIyclhyZIlhIaGOvGdXHq01jz00EPs3LmTBx98kNTUVI8cWzabRkZGaGtro7OzE3AsZxAQEEBUVBQ33XQTra2tRrfoc5Vr/f39SU9PZ/ny5eTm5tLY2EhbW5uRs+rr6wkMDCQ0NFTKyF8hPDycG2+8cVqjeFdXF4cOHcJqtVJQUEBTUxN5eXl89NFHPPbYYwCEhoYSGRlJTEwMq1evZtmyZVx55ZUkJCTM2ZIUHlcBHBoaYsOGDTQ1NfH666+TmJjo7JCEOG9eXl4kJiaSmJjI9ddfDzhaj/bt20dNTQ1NTU10d3czNjbGrl27aGtro6ysjImJCXbs2GHsx8/Pj9DQUAICAggLC2PBggUsWLCANWvWYLPZCA0NJSgo6IILSfKFP/OsViu33HILJpOJ119/Xb6MXZDVaqWpqYnOzk4OHTpEfX09fn5+XHbZZZhMJjZv3nze14aXlxcxMTHExMSQmZnJrl27jM98aGjI2L/FYiEsLIyQkBCXXaNrtplMJhYsWMCGDRvYsGED4MiHn332GYWFhZSXl1NdXU1nZye5ubnGOmDguFs4f/58oqOjiYuLIzExkbGxMcbHxz2ycv11jI6O8rOf/YznnnuOnJwcvvWtbzk7pEuGzWbjyJEjdHZ2GndUly1bZvQCmZrdMjY2lvb29gvev7e3t1GmmJiYABx3cktLS7HZbHh7ezN//nzMZjP9/f0EBQXN3Ju7RIWHh7Nx48ZplUK73U5LSwsAb7zxBlVVVbS2tlJdXc0nn3xi/K2XlxeRkZEEBwcTGRlJZGQk0dHRWK1Wo6vuTI199qgKYGtrKzfddBMlJSX89re/ZdGiRTIrorgkWCwWFi9ebLTenZ54JiYmaGpqwsfHhwMHDtDX10dbWxstLS3U19dTXl7O/v370VpPqySCo4AVHByMn58fwcHBREVFsXDhQlJTU2loaCA4OJiQkBCCg4MZGxuju7v7C4tRi4tz6NAh+vr6ePDBB2lububVV19lyZIlzg7Lo2mtGRoaoqenB6vVaozh7erqorq6mrCwMK644gqSkpJmpOV8qpv3FD8/P2677TZ2795txNDe3s4111wDOCYtOHToEIGBgYyPj9PQ0MC8efMICwv7WnG4m4CAAGM5CXDkw8LCQhobG2lqajK6vVVUVHDkyBEKCgoAeOKJJwDH7ImBgYHMnz+f9PR0p70PVzX13TIyMsKePXt47bXXKC8vZ+PGjWzdupWrrrrKyRG6H601fX19NDQ0YLVajQYks9lMc3MzAQEBRg+C9evXz0pj69SdJ5PJxKZNm2hvb6erq4uuri6KioqMnkWjo6Pk5eXR3t5OUFAQAQEBMizhHEwmk7Ee9OnjAu12O01NTTQ1NdHR0UFjYyPV1dVUVVVRVVXF6OgoAE899RTgqCAGBgYSFhb2tW9gzWoFUCm1CXgGMAM7tNZPnPG6D/AXIAvoBr6jta6b6Ti01vzzn//khz/8IR0dHfz+97/nF7/4hdytEB7BbDaTkJBAVlbWtBnZTq8kjo2N0dHRQXR0NEVFRXR2dlJTU0N3dzdaaxoaGujv7+fo0aPs37+foaGhr/yf3t7e+Pj4GA9fX198fX257LLL0Frj5eVltLSnpKTM6vs/G1fJTV9l165dPP3003R1dfHYY49x1113zeW/92h2u53h4WG8vLywWCz09vaSm5vL4cOHjZnzzGaz0SKfmJhIUlISXl5eZGZmztp3i8lkIjw8fNp4YZvNRmJiIiUlJSilCAwMZGBggLKyMmpqagCMxdmbmpo4fvw4vr6++Pj4MD4+Tm1tLQkJCfj4+DAyMkJvby8jIyN4e3vPWVekueDt7U1ycjLJycnAf/Kf1hqr1crJkyexWCyUlpbS2tpKWVkZLS0tVFZWznmsc5GftNZGwbO5uZmKigqGhoYYHR017oT6+/vj6+vL6OgoPT09aK2x2+20tbXR2tpKTU0NIyMjJCUlsXXrVtatW8f69euN7cQXaa05deqU0Y0T4PDhwzQ3NxMdHU19fT1ms9lYlkwpxaZNm+a854fFYiEuLs5YVD4tLY3i4mIATp06RXt7O2VlZcb2SinCw8MBGBwcpLGxEV9fX8LDw+nt7TVyjvRgmc5kMhmzIU85vWzW3d1NQ0MDQUFBFBcX09LSQkVFBZ2dnRw5cuRr/e9ZqwAqpczA/wM2AE3AAaVUrta6/LTNvg9YtdapSqnvAk8C37nQ/3W2ySg+++wzysvLKSgoID8/n/r6etLS0nj55ZfZtGnTBe1LiEudxWIxWqfO7F5wtglfli9fzscff4zVaqW3txer1UpQUBANDQ2Mjo7S2dlJS0sLp06d4tSpU4yMjNDX18fRo0fp7+9ncHCQU6dOzeVbNDg7N31ZjhkdHaWkpISdO3fy9ttvU1FRQVhYGM8//zwrV678wv4kN50/rTWjo6NMTEzQ1dVFR0cHY2NjxljYkZERSktLjfP18OHD1NXVkZGRQVJSklERio2NJTQ0lNDQUAIDA4mMjKS5udmpSzZ4eXkZ8YWGhrJ69WrAsdTM3r17GR4eNrpteXt7ExAQwMjICIODg0YXyakW6bq6OsrKyqivrwcclVxvb2+WLVtmvN7a2orZbGZwcJDq6mrMZjNLly5FKUVXVxdDQ0OYTCajQNje3s7ChQsBR8Gwq6uL/v5+Y5vTK5lTE1yBo+eC3W5HKTWrhUalFGFhYVx11VVnvT4zMzPndLmJ2cpPAwMDPP7445SUlFBZWUl1dfU5G/LMZrPxsNvtmEwm41gEBwezbt06br31VrKysvjTn/5kTIIBjrUvo6KiiI6ONrpGV1dXU1dXZ+wzPT2dkJAQenp6jN9PTExw4sQJY5kQX19f+vv7aW5uNv5OKUVbWxs2mw0vLy9GRkYYHR3FarUyMDBgnDNTd6VsNptxLo2PjzMxMYFSasY+V5vNZlSax8bGqK2tpa6uzrhDc+zYMT777DPju9But1NbW8vSpUsBxzjV2NhYsrOzaW5uJjAwcFpsrlBp8vPzw8fHB3B0O73++uvZu3cvAwMDDA0NMTQ0RGhoKD09PfT19RnnQU9Pj5FP1q1bR3h4OLW1tXzyySdYLBYsFgve3t6MjY1ht9vx9fWdVqZoaWmZdg6C47tyfHwcpRRDQ0PG2qkWi8U5B2cWTa2fmpWV9YXG8qysrK91bszmHcBVQJXWugZAKfU3YAtwehLbAjw6+fwd4FmllNJfMYVheXk5S5YswW63G18Q4+PjWCwW7Ha7cTIMDAwYfxMfH8+OHTtYtmwZXl5eHDp0SApPQnwN3t7eLFy40CjYwdkrN2c6s0/80NAQ8fHxc30XcFZyE8CRI0eIi4tDa43WmrGxMUJDQ40C89QXF2BsEx4eTn9/v3G3FRxfsNdeey033HCDUfguKytj9+7dxv+a6iqSk5ODUorKykoaGhqorq42tjGZTEauKy8vp62tDXAM8K+ursZisbB27VrAMVPl559/Pu39dHZ2EhERAThaqWtra6ftPzAw0Nh/UVGRsf5UTU0NNTU1xvgUgA8//JCSkpJp+x8eHjYKaXv27GF0dHTaDLZTi+8C5OXlGcfu5MmT1NTUEBMTYxyfN998k+rqauO7YWJigvHxcby9vZmYmGDnzp0AHD9+3CiQTHWZNplM9PX14efnR0REBBkZGYSFhRkt2oGBgVx77bVu1WtEKWXceZ9ytmu2sLDQWFMyJiaG2NhY407n1GOqgmuz2RgeHja6lTc2NjIxMWF8Bg0NDdTV1Rn7b29vp7m5mVtvvRVwnIOVlZXG8QdH4XfNmjWA4xyaOkdLS0upr68nICCAG2+8EYAPPviAwsJCo4B/9OhRent7yc7OBqCgoIDKykpqa2uN/U8VngDy8/MZHh421kerqakhIiLCaGB5//33KSsrMwpVSiln3Mmalfx04sQJtm3bZnTrj4qKIjk5mZtvvpnAwEDee+89YmNj6e7uNirnd9xxB3FxcQwMDPDiiy9it9sJCQmhtbUVu93O2rVrWb16NYmJicTHxxvlMbvdTn9/PwEBAYCjh0ljYyNjY2M0Nzdjt9ux2+0sWLCAkJAQ2tvbycvLw263c+DAATo6OgC48847WbRoESUlJfz5z3+e9n4iIiLYtGkT8fHxFBcXk5eXR0RExLQ7a9///vcBjBsB4Bif1dXVBcBPfvITQkJCeOutt3jjjTemFaYjIiK4++678fPzY+fOnZSUlBgNGlprTCYTDz30EADvvfcex44dM/4+PDycvr4+HnroIUwmE4WFhZSUlODt7Y3FYsHX13daBW94eBir1Uppaalx1z4oKMjoSltYWMjg4KCxfU1NDVar1Tiv9+7d+4XcPDo6alTYPv/8c8bGxqirqzO2iYyMNCqgu3fvnnae19fXc+rUKWNG8d27d3/hu2Wqd4/ZbObYsWPGfqa2WblyJVFRUcTFxbF9+3bGx8cpLi7G29ubkJAQhoeH8fb2pru7m8rKSo4dO0ZkZCS+vr60tLQwMjJCT0+P8b0G0N/fT2Njo5F3AgMDOXjwIPv27QNg1apVmM1mdu3aRWtrK4GBgQwODtLS0sLChQtJSEgAHA1Z/v7+VFRUGHf8zWYzK1asoLm5mU8++YSwsDB8fX0pLCykqKgIi8ViDMH45JNPjHrHVINicnKycb5PNarV19ezb98+urq68PPzIzU1lebmZj788EMSExMxm820tbXR3d3NvHnzjN4JO3fuNL4XW1pasFqtLF261DhnTpw4waJFi2hububrUDM9XbyxY6XuADZprR+Y/Plu4Cqt9U9P2+bY5DZNkz9XT27Tdca+fgD8YPLHy4G575NxbuFA1zm3ci53iBEkzpnmDnFerrWek6lJZzI3Tb7m6vnJHT5/kDhnkjvECO4Tp1vmJzfITeAe54A7xAgS50xyhxjha+Qmt5gERmu9Hdju7Di+ilLqoNbapRcbcocYQeKcae4Qp1LqoLNjuFiunp/c4fMHiXMmuUOM4F5xOjuGi+HquQnc4xxwhxhB4pxJ7hAjfL3cNJsd25uBuNN+jp383Vm3UUp5AcE4BjQLIcRskdwkhHBVkp+EELNuNiuAB4BFSqkkpZQF+C6Qe8Y2ucC9k8/vAPLONcZGCCG+JslNQghXJflJCDHrZq0LqNbappT6KfARjqmMX9Zalyml/hs4qLXOBV4CXlVKVQE9OBKdu3LpbhaT3CFGkDhnmjvEOWcxSm5yWRLnzHGHGEHi/ALJTy7JHWIEiXMmuUOM8DXinLVJYIQQQgghhBBCuJa5W9xGCCGEEEIIIYRTSQVQCCGEEEIIITyEVAAvgFJqk1KqUilVpZT6xVlev08p1amUKpl8POCkOF9WSnVMrhV0tteVUuqPk+/jiFIq0wVjvFYp1Xfasfz1XMc4GUecUupTpVS5UqpMKfW/z7KNU4/necbo9OOplPJVShUppUon4/y/Z9nGRyn15uSx3K+USpzrON2VO+Qnd8hNk3G4fH5yh9x0AXG6wvGU/DRLJDfNHHfITZNxuHx+8vjcpLWWx3k8cAzGrgaSAQtQCqSdsc19wLMuEOs6IBM49iWv3wT8C1DAamC/C8Z4LfC+CxzLKCBz8nkgcOIsn7tTj+d5xuj04zl5fAImn3sD+4HVZ2zzY+CFyeffBd509jngDg93yU/ukJvOM05XuJ5cPjddQJyucDwlP83OcZXcNLdxOv1amozD5fOTp+cmuQN4/lYBVVrrGq31GPA3YIuTYzorrfUeHDODfZktwF+0QyEQopSKmpvoHM4jRpegtW7VWh+efD4AVAAxZ2zm1ON5njE63eTxGZz80XvyceYsVFuAVyafvwNcr5RScxSiO3OL/OQOuQncIz+5Q266gDidTvLTrJHcNIPcITeBe+QnT89NUgE8fzFA42k/N3H2E+X2yVvZ7yil4s7yuis43/fibFdP3vL+l1JqqbODmbylnoGj9eV0LnM8vyJGcIHjjkEKSQAABIRJREFUqZQyK6VKgA7gY631lx5LrbUN6APmz22UbulSyU8ucy2dB6dfT1PcITeB5CcPJblp7jn9WjqdO+QnT8xNUgGcWf8EErXWy4GP+U9tXFy4w0CC1noF8CfgXWcGo5QKAP4O/B+tdb8zY/ky54jRJY6n1npCa70SiAVWKaWWOSMODyX5aea4xPUE7pGbQPKT+EqSm2aOS1xLU9whP3lqbpIK4PlrBk5vlYqd/J1Ba92ttR6d/HEHkDVHsV2oc74XZ9Na90/d8tZa7wS8lVLhzohFKeWNIzn8VWv9j7Ns4vTjea4YXel4TsbQC3wKbDrjJeNYKqW8gGCge26jc0uXSn5y+rV0PlzlenKH3ASSnzyc5KY55ErXkjvkJ0/OTVIBPH8HgEVKqSSllAXHIMvc0zc4o+/yrTj6E7uiXOCeyRmYVgN9WutWZwd1OqVU5FT/ZaXUKhzn6px/0U7G8BJQobX+w5ds5tTjeT4xusLxVEpFKKVCJp/7ARuA42dslgvcO/n8DiBPa31mX3fxRZdKfnL53AQucz25fG4CyU9CctNccoVrafJ/u3x+8vTc5DXTgV6qtNY2pdRPgY9wzGr1sta6TCn138BBrXUu8L+UUrcCNhyDdO9zRqxKqTdwzFwUrpRqAn6DY9AoWusXgJ04Zl+qAoaB+10wxjuAHymlbMAp4LtO+qJdA9wNHJ3sfw3wSyD+tFidfTzPJ0ZXOJ5RwCtKKTOOJPqW1vr9M66hl4BXlVJVOK6h785xjG7JXfKTO+Sm84zTFa4nd8hN5xunKxxPyU+zQHLTnMfpCtcSuEd+8ujcpKTxSgghhBBCCCE8g3QBFUIIIYQQQggPIRVAIYQQQgghhPAQUgEUQgghhBBCCA8hFUAhhBBCCCGE8BBSARRCCCGEEEIIDyEVQCGEEEIIIYTwEFIBFEIIIYQQQggPIRVA4TKUUv+llHph8rm3UupVpdQrSilvZ8cmhPBckpuEEK5K8pO4GFIBFK4kHTiilAoC/gU0aK3v1VqPOzkuIYRnk9wkhHBVkp/EBZMKoHAly4EeYDfwttZ6m3PDEUIIQHKTEMJ1SX4SF0xprZ0dgxAAKKWsgA34ntb6n86ORwghQHKTEMJ1SX4SF8PL2QEIAaCUigMGgZNA1Gm/XwpsA7qAJq31/zgnQiGEJ5LcJIRwVZKfxMWSCqBwFelAKbAVKFRKHdBaFwM3Aq9qrf/l1OiEEJ5KcpMQwlVJfhIXRcYAClexHDiqtW4FHgDeVEoFAy8Ba5VSLyml/supEQohPJHkJiGEq5L8JC6KVACFq0gHjgJorT8G3gJe1lr3a623aa2/D9yqlJJzVggxlyQ3CSFcleQncVFkEhjh0pRS3wQ24hjgPKK1/rmTQxJCCMlNQgiXJflJnItUAIUQQgghhBDCQ8gtYSGEEEIIIYTwEFIBFEIIIYQQQggPIRVAIYQQQgghhPAQUgEUQgghhBBCCA8hFUAhhBBCCCGE8BBSARRCCCGEEEIIDyEVQCGEEEIIIYTwEFIBFEIIIYQQQggP8f8BuM6PzAUIkGYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x216 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,3, figsize=(15,3))\n",
    "for i, m in enumerate(models):\n",
    "    plot_mixture(m, data=X, ax=ax[i], l=0.05, u=3, bins=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems more parsimonious to use the two-component model, so I will use that one. I will  now select components from the first component."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>AlignmentCoverage</th>\n",
       "      <th>AlignmentIdentity</th>\n",
       "      <th>AlignmentLength</th>\n",
       "      <th>AlignmentLengthStripped</th>\n",
       "      <th>Distance</th>\n",
       "      <th>Family</th>\n",
       "      <th>Ka</th>\n",
       "      <th>Ks</th>\n",
       "      <th>Node</th>\n",
       "      <th>Omega</th>\n",
       "      <th>Outlier</th>\n",
       "      <th>Paralog1</th>\n",
       "      <th>Paralog2</th>\n",
       "      <th>WeightOutliersExcluded</th>\n",
       "      <th>WeightOutliersIncluded</th>\n",
       "      <th>log(Ks)</th>\n",
       "      <th>p_component1</th>\n",
       "      <th>p_component2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>AT4G25420__AT5G51810</th>\n",
       "      <td>0.74008</td>\n",
       "      <td>0.78731</td>\n",
       "      <td>1512.0</td>\n",
       "      <td>1119.0</td>\n",
       "      <td>0.29210</td>\n",
       "      <td>GF_000197</td>\n",
       "      <td>0.1211</td>\n",
       "      <td>0.7495</td>\n",
       "      <td>13.0</td>\n",
       "      <td>0.1615</td>\n",
       "      <td>False</td>\n",
       "      <td>AT4G25420</td>\n",
       "      <td>AT5G51810</td>\n",
       "      <td>1.00</td>\n",
       "      <td>1.00</td>\n",
       "      <td>-0.288349</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>5.872964e-09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AT1G15540__AT1G80320</th>\n",
       "      <td>0.62698</td>\n",
       "      <td>0.59599</td>\n",
       "      <td>1512.0</td>\n",
       "      <td>948.0</td>\n",
       "      <td>0.73190</td>\n",
       "      <td>GF_000197</td>\n",
       "      <td>0.3517</td>\n",
       "      <td>68.8938</td>\n",
       "      <td>20.0</td>\n",
       "      <td>0.0051</td>\n",
       "      <td>True</td>\n",
       "      <td>AT1G80320</td>\n",
       "      <td>AT1G15540</td>\n",
       "      <td>0.00</td>\n",
       "      <td>1.00</td>\n",
       "      <td>4.232566</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>4.218759e-42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AT1G54280__AT3G13900</th>\n",
       "      <td>0.91834</td>\n",
       "      <td>0.84937</td>\n",
       "      <td>4041.0</td>\n",
       "      <td>3711.0</td>\n",
       "      <td>0.12266</td>\n",
       "      <td>GF_000206</td>\n",
       "      <td>0.0581</td>\n",
       "      <td>1.0479</td>\n",
       "      <td>13.0</td>\n",
       "      <td>0.0554</td>\n",
       "      <td>False</td>\n",
       "      <td>AT3G13900</td>\n",
       "      <td>AT1G54280</td>\n",
       "      <td>1.00</td>\n",
       "      <td>1.00</td>\n",
       "      <td>0.046788</td>\n",
       "      <td>0.999935</td>\n",
       "      <td>6.534271e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AT1G17500__AT1G72700</th>\n",
       "      <td>0.90052</td>\n",
       "      <td>0.84776</td>\n",
       "      <td>4041.0</td>\n",
       "      <td>3639.0</td>\n",
       "      <td>0.12898</td>\n",
       "      <td>GF_000206</td>\n",
       "      <td>0.0254</td>\n",
       "      <td>0.9952</td>\n",
       "      <td>14.0</td>\n",
       "      <td>0.0256</td>\n",
       "      <td>False</td>\n",
       "      <td>AT1G72700</td>\n",
       "      <td>AT1G17500</td>\n",
       "      <td>1.00</td>\n",
       "      <td>1.00</td>\n",
       "      <td>-0.004812</td>\n",
       "      <td>0.999982</td>\n",
       "      <td>1.819869e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AT1G72700__AT3G13900</th>\n",
       "      <td>0.91017</td>\n",
       "      <td>0.74008</td>\n",
       "      <td>4041.0</td>\n",
       "      <td>3678.0</td>\n",
       "      <td>0.32305</td>\n",
       "      <td>GF_000206</td>\n",
       "      <td>0.0796</td>\n",
       "      <td>2.4245</td>\n",
       "      <td>15.0</td>\n",
       "      <td>0.0328</td>\n",
       "      <td>False</td>\n",
       "      <td>AT3G13900</td>\n",
       "      <td>AT1G72700</td>\n",
       "      <td>0.25</td>\n",
       "      <td>0.25</td>\n",
       "      <td>0.885625</td>\n",
       "      <td>0.040212</td>\n",
       "      <td>9.597877e-01</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                      AlignmentCoverage  AlignmentIdentity  AlignmentLength      ...        log(Ks)  p_component1  p_component2\n",
       "AT4G25420__AT5G51810            0.74008            0.78731           1512.0      ...      -0.288349      1.000000  5.872964e-09\n",
       "AT1G15540__AT1G80320            0.62698            0.59599           1512.0      ...       4.232566      1.000000  4.218759e-42\n",
       "AT1G54280__AT3G13900            0.91834            0.84937           4041.0      ...       0.046788      0.999935  6.534271e-05\n",
       "AT1G17500__AT1G72700            0.90052            0.84776           4041.0      ...      -0.004812      0.999982  1.819869e-05\n",
       "AT1G72700__AT3G13900            0.91017            0.74008           4041.0      ...       0.885625      0.040212  9.597877e-01\n",
       "\n",
       "[5 rows x 18 columns]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = models[1]\n",
    "anchors_ks = anchors_ks[anchors_ks[\"Ks\"] > 0]  # filter out Ks = 0, as log is -inf!\n",
    "ks_anchors_mix = get_component_probabilities(df=anchors_ks, model=model)\n",
    "ks_anchors_mix.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4332 paralogs were selected\n"
     ]
    }
   ],
   "source": [
    "paralogs = set(ks_anchors_mix[ks_anchors_mix['p_component1'] > 0.95]['Paralog1'])\n",
    "paralogs = paralogs | set(ks_anchors_mix[ks_anchors_mix['p_component1'] > 0.95]['Paralog2'])\n",
    "print('{} paralogs were selected'.format(len(paralogs)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. GO enrichment analysis\n",
    "\n",
    "I have the genes of interest, _i.e._ the ones I hope to resemble paralogs retained from a WGD represented by the peak I modeled with the Bayesian Gaussian mixture model. Now I will perform GO enrichment on the selected paralogs using the `goenrich` library. To this end I use the whole gene ontology obo graph (retrieved at http://purl.obolibrary.org/obo/go/go-basic.obo and stored in the directory `~/godb`) and the GO annotatoin for _A. thaliana_ obtained from the GO consortium website (http://geneontology.org/gene-associations/gene_association.tair.gz, also in the directory `~/godb`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Set up the GO enrichment analysis**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/arzwa/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py:2961: DtypeWarning: Columns (3) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  exec(code_obj, self.user_global_ns, self.user_ns)\n"
     ]
    }
   ],
   "source": [
    "ontology = goenrich.obo.ontology('/home/arzwa/godb/go-basic.obo')\n",
    "gene2go = goenrich.read.sgd('/home/arzwa/godb/gene_association.tair.gz')\n",
    "values = {k: set(v) for k,v in gene2go.groupby('go_id')['db_object_name']}\n",
    "background_attribute = 'sgd'\n",
    "goenrich.enrich.propagate(ontology, values, background_attribute)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Perform enrichment analysis and get significant results**\n",
    "\n",
    "The underlying test is the `hypergeometric test hypergeom.sf(x, M, n, N)`. All $p$-values are corrected for multiple testing using the Benjamini-Hochberge method before a significance cut-off of $\\alpha=0.05$ is applied."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = goenrich.enrich.analyze(ontology, paralogs, background_attribute)\n",
    "df = df.dropna()\n",
    "enriched = df[df['rejected'] == 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 327 enriched BP terms\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>term</th>\n",
       "      <th>q</th>\n",
       "      <th>rejected</th>\n",
       "      <th>p</th>\n",
       "      <th>x</th>\n",
       "      <th>n</th>\n",
       "      <th>M</th>\n",
       "      <th>N</th>\n",
       "      <th>name</th>\n",
       "      <th>namespace</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>24535</th>\n",
       "      <td>GO:0048468</td>\n",
       "      <td>8.966824e-10</td>\n",
       "      <td>1.0</td>\n",
       "      <td>6.625732e-13</td>\n",
       "      <td>80</td>\n",
       "      <td>320</td>\n",
       "      <td>39891</td>\n",
       "      <td>4332</td>\n",
       "      <td>cell development</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4131</th>\n",
       "      <td>GO:0005976</td>\n",
       "      <td>8.966824e-10</td>\n",
       "      <td>1.0</td>\n",
       "      <td>5.660499e-13</td>\n",
       "      <td>104</td>\n",
       "      <td>465</td>\n",
       "      <td>39891</td>\n",
       "      <td>4332</td>\n",
       "      <td>polysaccharide metabolic process</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31150</th>\n",
       "      <td>GO:0071396</td>\n",
       "      <td>7.478127e-09</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.209516e-12</td>\n",
       "      <td>90</td>\n",
       "      <td>397</td>\n",
       "      <td>39891</td>\n",
       "      <td>4332</td>\n",
       "      <td>cellular response to lipid</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33743</th>\n",
       "      <td>GO:0090627</td>\n",
       "      <td>3.298399e-08</td>\n",
       "      <td>1.0</td>\n",
       "      <td>6.499309e-11</td>\n",
       "      <td>48</td>\n",
       "      <td>163</td>\n",
       "      <td>39891</td>\n",
       "      <td>4332</td>\n",
       "      <td>plant epidermal cell differentiation</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33678</th>\n",
       "      <td>GO:0090558</td>\n",
       "      <td>4.277078e-08</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.481208e-11</td>\n",
       "      <td>67</td>\n",
       "      <td>272</td>\n",
       "      <td>39891</td>\n",
       "      <td>4332</td>\n",
       "      <td>plant epidermis development</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9059</th>\n",
       "      <td>GO:0016051</td>\n",
       "      <td>6.029326e-08</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.485056e-10</td>\n",
       "      <td>84</td>\n",
       "      <td>379</td>\n",
       "      <td>39891</td>\n",
       "      <td>4332</td>\n",
       "      <td>carbohydrate biosynthetic process</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4816</th>\n",
       "      <td>GO:0006812</td>\n",
       "      <td>2.998255e-07</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.204314e-10</td>\n",
       "      <td>89</td>\n",
       "      <td>425</td>\n",
       "      <td>39891</td>\n",
       "      <td>4332</td>\n",
       "      <td>cation transport</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19256</th>\n",
       "      <td>GO:0042546</td>\n",
       "      <td>3.751221e-07</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.293524e-09</td>\n",
       "      <td>58</td>\n",
       "      <td>234</td>\n",
       "      <td>39891</td>\n",
       "      <td>4332</td>\n",
       "      <td>cell wall biogenesis</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>526</th>\n",
       "      <td>GO:0000904</td>\n",
       "      <td>4.628327e-07</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.823971e-09</td>\n",
       "      <td>58</td>\n",
       "      <td>236</td>\n",
       "      <td>39891</td>\n",
       "      <td>4332</td>\n",
       "      <td>cell morphogenesis involved in differentiation</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>189</th>\n",
       "      <td>GO:0000271</td>\n",
       "      <td>5.615480e-07</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.766246e-09</td>\n",
       "      <td>53</td>\n",
       "      <td>209</td>\n",
       "      <td>39891</td>\n",
       "      <td>4332</td>\n",
       "      <td>polysaccharide biosynthetic process</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21821</th>\n",
       "      <td>GO:0045488</td>\n",
       "      <td>8.095366e-07</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4.386651e-09</td>\n",
       "      <td>47</td>\n",
       "      <td>177</td>\n",
       "      <td>39891</td>\n",
       "      <td>4332</td>\n",
       "      <td>pectin metabolic process</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9057</th>\n",
       "      <td>GO:0016049</td>\n",
       "      <td>8.344577e-07</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4.932755e-09</td>\n",
       "      <td>80</td>\n",
       "      <td>380</td>\n",
       "      <td>39891</td>\n",
       "      <td>4332</td>\n",
       "      <td>cell growth</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7445</th>\n",
       "      <td>GO:0010393</td>\n",
       "      <td>8.666898e-07</td>\n",
       "      <td>1.0</td>\n",
       "      <td>5.336760e-09</td>\n",
       "      <td>47</td>\n",
       "      <td>178</td>\n",
       "      <td>39891</td>\n",
       "      <td>4332</td>\n",
       "      <td>galacturonan metabolic process</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31420</th>\n",
       "      <td>GO:0071669</td>\n",
       "      <td>1.240750e-06</td>\n",
       "      <td>1.0</td>\n",
       "      <td>8.207174e-09</td>\n",
       "      <td>64</td>\n",
       "      <td>282</td>\n",
       "      <td>39891</td>\n",
       "      <td>4332</td>\n",
       "      <td>plant-type cell wall organization or biogenesis</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24905</th>\n",
       "      <td>GO:0048868</td>\n",
       "      <td>1.240750e-06</td>\n",
       "      <td>1.0</td>\n",
       "      <td>8.251291e-09</td>\n",
       "      <td>48</td>\n",
       "      <td>186</td>\n",
       "      <td>39891</td>\n",
       "      <td>4332</td>\n",
       "      <td>pollen tube development</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             term             q         ...                                                     name           namespace\n",
       "24535  GO:0048468  8.966824e-10         ...                                         cell development  biological_process\n",
       "4131   GO:0005976  8.966824e-10         ...                         polysaccharide metabolic process  biological_process\n",
       "31150  GO:0071396  7.478127e-09         ...                               cellular response to lipid  biological_process\n",
       "33743  GO:0090627  3.298399e-08         ...                     plant epidermal cell differentiation  biological_process\n",
       "33678  GO:0090558  4.277078e-08         ...                              plant epidermis development  biological_process\n",
       "9059   GO:0016051  6.029326e-08         ...                        carbohydrate biosynthetic process  biological_process\n",
       "4816   GO:0006812  2.998255e-07         ...                                         cation transport  biological_process\n",
       "19256  GO:0042546  3.751221e-07         ...                                     cell wall biogenesis  biological_process\n",
       "526    GO:0000904  4.628327e-07         ...           cell morphogenesis involved in differentiation  biological_process\n",
       "189    GO:0000271  5.615480e-07         ...                      polysaccharide biosynthetic process  biological_process\n",
       "21821  GO:0045488  8.095366e-07         ...                                 pectin metabolic process  biological_process\n",
       "9057   GO:0016049  8.344577e-07         ...                                              cell growth  biological_process\n",
       "7445   GO:0010393  8.666898e-07         ...                           galacturonan metabolic process  biological_process\n",
       "31420  GO:0071669  1.240750e-06         ...          plant-type cell wall organization or biogenesis  biological_process\n",
       "24905  GO:0048868  1.240750e-06         ...                                  pollen tube development  biological_process\n",
       "\n",
       "[15 rows x 10 columns]"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BP = enriched[enriched['namespace'] == 'biological_process']\n",
    "BP = BP.sort_values('q', ascending=True)\n",
    "print(\"There are {} enriched BP terms\".format(len(BP)))\n",
    "BP.head(n=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$M$ is the total number of genes (background set), $N$ the number of genes in the ste of interest, $n$ the number of genes annotated with the particular term in the background set and $x$ the number of genes with the term of interest annotated in the test set of interest.\n",
    "\n",
    "Notice how many development and defense related terms appear."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
